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SECTION  I 

INTRODUCTION 


Vibration  issues  and  the  associated  fatigue  accumulation  often  play  a  critical  role 
in  the  design  of  structural  components.  This  observation  is,  in  particular,  true  in  the 
context  of  future  supersonic/hypersonic  vehicles  such  as  the  planned  National  Aerospace 
Space  Plane  (NASP)  the  skin  of  which  will  be  subjected  to  especially  harsh  operating 
conditions,  e.g.  surface  temperatures  possibly  exceeding  3000°F,  severe  acoustic  loading 
from  the  engine  exhaust,  etc.  The  design  of  the  surface  panels  appears  in  this  light  quite 
challenging  especially  since  it  is  not  only  required  to  consider  each  of  these 
environmental  factors  by  itself  but  it  is  also  necessary  to  account  for  their  combined 
effects.  As  an  example  of  loading  interaction,  consider  the  response  of  panels  to  the 
thermal  and  acoustic  excitations  mentioned  above  and  note  first  that  the  increase  of 
surface  temperature  will  produce  compressive  stresses  in  the  panels  since  the  extension 
that  they  would  naturally  undergo  is  prevented  by  their  supports.  Further,  the  magnitude 
of  these  stresses,  in  direct  relation  to  the  very  large  temperature  changes,  is  likely  to 
produce  the  buckling  of  some  panels. 

To  sustain  such  high  surface  temperatures,  specially  designed  structural  materials, 
such  as  ceramic  matrix  composites,  will  have  to  be  used.  The  refractory  nature  of  these 
materials  will  prevent  the  heat  conduction  through  the  panel  and  thus  will  lead  to  a  severe 
temperature  gradient  in  that  direction.  The  corresponding  compressive  stresses  will  exhibit 
a  similar  sharp  variation  through  the  thickness  which  will  result  not  only  in  a  normal  force 
but  also  in  a  moment  that  increases  the  likelihood  of  buckling  (see  Ng,  1988,  1989,  Lee, 
1993,  1997,  Vaicaitis,  1994,  Moorthy  et  al.  1995) 

Of  primary  importance  in  the  fatigue  damage  accumulation  process  is  the  cycling 
of  the  stresses  and,  thus,  buckling  represents  a  specially  acute  problem  if  the  panel 
oscillates  from  its  buckled  state  to  its  normal  configuration,  or  from  one  buckled  state  to 
another.  This  latter  mechanism,  often  referred  to  as  snap-through  or  oil-canning,  is 
especially  likely  when  the  panel  is  subjected,  in  addition  to  the  thermal  loading,  to  a 
random  transverse  excitation  such  as  the  acoustic  loading  from  the  engine  exhaust.  The 
corresponding  response  of  the  panel  consists  of  random  fluctuations  around  the  buckled 
position  but  also  often  includes  large  excursions  from  this  configuration,  so  large  in  fact 
that  the  panel  may  snap-through  to  the  other  buckled  state.  Accordingly,  the  goal  of  the 
present  investigation  is  the  formulation  and  preliminary  assessment  of  a  methodology  for 
the  prediction  of  the  fatigue  damage  accumulated  in  a  panel  due  to  both  fluctuations 
around  the  buckled  states  and  the  snap-throughs. 


1 


SECTION  2 

PANEL  STRUCTURAL  DYNAMIC  MODELING 

The  structural  dynamic  modeling  of  the  panel  can  be  decomposed  into  two  major 
parts:  the  composite  plate  modeling  and  the  derivation  of  a  simplified,  one-mode,  model 
both  of  which  are  discussed  below. 

2.1  Composite  Plate  Modeling 

Recent  efforts  in  the  structural  dynamics  of  composite  plates  (see  for  example 
Chattopadhyay  and  Gu,  1994)  have  emphasized  the  appropriateness  of  the  higher  order 
plate  theory  as  proposed  by  Reddy  (1987).  Specifically,  the  displacement  field  inside  the 
plate  are  selected  in  the  form  (see  Fig.  2.1) 

A  z 


Figure  2.1  Composite  plate  modeling 


u(x,  y,  z,  t)  =  u0(x,  y,t)-z 


dw0(x,  y,  t) 


dx 


+  z 


y,t) 


(i) 


v(x,  y,  z,  0  =  v0(x,  y,t)-z 


dy 


(2) 


w(x,  y,  z,  t )  =  w0  (x,  y,  t)  (3) 

where  «,  v,  and  w  are  the  displacements  along  the  x,  y,  and  z  axes,  respectively  and  h 
denotes  the  plate  thickness.  Before  proceeding  further,  it  is  important  to  assess  the  origin 
of  each  term  present  in  the  above  equations. 

First,  uq(x,  y,  t)  and  v0(x,  y,  t )  represent  the  time-varying  in-plane  displacements 

of  the  mid-plane.  These  displacements  are  produced  by  both  the  thermal  effects,  which 
include  an  overall  increase  of  the  panel  temperature  and  temperature  gradients  along  and 
across  the  plate,  and  the  mid-plane  stretch  associated  with  the  large  transverse 
displacements.  Note  that  both  of  these  effects  play  fundamental  roles  in  the  present 
buckling/postbuckling  analysis  and  thus  must  absolutely  be  included  in  the  formulation  of 
the  plate  equations.  Clearly,  the  thermal  effects  induce  compressive  in-plane  stresses  that 
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/ 


are  reflected  by  an  apparent  softening  of  the  plate  in  the  transverse  direction  and 
consequently  facilitate  the  plate  buckling.  Once  buckling  is  initiated,  the  transverse 
displacements  rapidly  increase  leading  to  a  stretching  of  the  midplate  which  restores  the 
stiffness  of  the  panel  and  allows  the  existence  of  a  stable  buckled  equilibrium  position. 
Clearly  then,  one  can  formally  write 

wo  =uo{T’wo)  and  vo=vo(^wo)  (4) 


where  the  dependence  of  mq  and  vo  on  w0  ls  intrinsically  nonlinear  as  it  describes  the 


mid-plane  stretching,  a  nonlinear  effect  of  the  transverse  motions. 

The  next  group  of  terms  present  in  the  in-plane  displacements  consists  of 

z  dwo(x’  0  an(j  z  d^o  (*>  Zli)  which  readily  are  recognized  as  the  first  order  plate 
dx  dy 

bending  components  and  are  linear  in  the  plate  transverse  coordinate.  Completing  the 


formulation  of 


u(x,y,t )  and  v(x,y,t)  are  the  terms 


and  z 


l-i  [L 
3  U 


^2 


cpv(jc,  y,  t)  which  represent  the  higher  order  plate  bending 


corrections.  Generally  speaking,  these  components  take  the  form  g(z)cp(x,  y,t)  and  are 
designed  to  account  for  a  nonlinear  distribution  of  the  displacements  across  the  thickness. 


The  selection  of  the  function  H(z),  H(z)  =  z 


i-if- 

jU, 


,  and  the  physical 


interpretation  of  the  functions  (p„  (x,  yj)  andcpv(x,  y,  t)  will  be  discussed  in  the  next 
section  in  connection  with  the  strains.  Clearly,  these  terms  are  primarily  produced  by  the 
transverse  motions  and  thus 

<P«=(Pw(wo)  and  <Pve<Pv(wo)  •  (5) 

It  will  be  shown,  under  the  assumption  of  a  symmetric  composite  layering,  that  cpw  and 


cpv  exhibit  a  linear  functional  dependence  with  respect  to  w0 . 

Considering  finally  the  transverse  displacement,  it  is  seen  that  the  variations  of 
this  quantity  across  the  thickness  have  been  neglected  in  accordance  with  (reasonably) 
thin  plate  assumptions.  Then,  w0(x,  y,  t)  denotes  the  time-dependent  transverse 
displacement  of  the  mid-plane,  produced  by  external  loading  and  influenced  by  the  in¬ 
plane  (membrane)  stresses  as  described  earlier. 

The  next  step  in  the  structural  dynamic  modeling  of  the  panel  is  the  introduction 
of  the  strains-displacements  relationships.  The  presence  of  an  in-plane  displacement  field 

of  order  (w0 )° ,  i.e.  the  one  corresponding  to  the  thermal  effects  only,  requires  the 

consideration,  even  for  the  small  amplitude  transverse  vibration  problem,  of  nonlinear 
terms  in  the  definition  of  the  strains.  In  the  present  analysis,  it  will  be  assumed  that  the  in¬ 
plane  displacements  are  substantially  smaller  than  their  transverse  counterparts  so  that  the 
second  order  terms  in  u  and  v  are  neglected  but  those  in  w  are  retained.  This  approach 
yields  the  von  Karman  strains 
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_du  }_ 
dx  2 


'dw'2 


dx 


dv  1 

y  a-  2 


^  dxv^2 


dy 


dy 


Sz  = 


dw 

dz 


(6) 


du 

dv 

f  dw' 

f  dw' 

1 

du 

dw 

1 

"dv 

dw 

_dy 

+  — 
dx 

+ 

UxJ 

UtJ 

Yxz- 2 

dz 

+  — 
dx _ 

Jyz  -  2 

dz 

+~dy. 

lxy~2 

i_  — • 

Introducing  the  displacement  field  given  by  Eq.  (l)-(3)  in  Eq.  (6)  yields  the  expressions 
of  the  strains  in  terms  of  the  five  basic  unknowns,  i.e.  u0(x,  y,  t),  v0(x,  y,t),  w0(x,  y,  t), 
<pu(x,y,t),  and  <pv(x,  y,t).  In  this  respect,  note  at  the  contrary  of  first  order  bending 
theory,  that  the  shears  strains  yxz  and  yyz  do  not  vanish  inside  the  plate.  In  fact,  one  has 
sz\ 2]  1 [  /r\2 

1  -  4 i  A J  <Pu(x,y’t)  and  lyz  =  2  1"4U 


Yxz  = 


9  v(x,y,t)  (7) 


from  which  it  is  then  seen  that  the  function  g(z )  was  selected  so  that  the  shear  strains  yxz 


and  yyz  vanish  at  both  the  top  and  the  bottom  of  the  plate  (z  =  ±  -  )  as  required  by  the  no 
shear  stress  boundary  conditions.  Further,  Eq.  (7)  implies  that  cpM(x,  y,i)  and  cpv(x,  y,t) 


are  in  fact  directly  related  to  the  mid-plane  (z  =  0)  shears  yxz  and  yyz .  The  higher  order 


plate  bending  displacement  field  given  by  Eq.  (l)-(3)  thus  includes  shear  effects. 

Turning  now  to  the  definition  of  the  stresses,  it  will  be  assumed  that  the  behavior 
of  the  panel  remains  linearly  elastic  during  its  entire  fatigue  life.  Modeling  further  each 
layer  of  the  composite  as  an  orthotropic  material  leads  to  the  stress-strain  relationships 


<jx< 

1 

Gy 

Sy  —  Cty  T 

^z 

>  =  Q. 

*z 

V 

Tx'z 

£x'z 

Ty'z 

Sy'z 

rx'y' 

£x'y' 

-Qa'T 


(8) 


where  (x\  /,  z)  denotes  the  frame  of  reference  for  each  layer  with  x'  aligned  with  the 
fibers.  Further,  ax>  and  ay  represent  the  coefficients  of  thermal  expansion  along  and 

across  fibers  and  T  =  T(x,  y,  z)  is  the  local  temperature.  Finally,  sxy  =  2  yxy , 
ex’z  =  2  yx'z  »  and  similarly  syz  =  2  yyz ,  and  Q  denotes  the  symmetric  elastic  constant 
matrix  of  the  orthotropic  layer,  i.e., 
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(9) 


Qw  Qm  0i3  ooo 

012  022  023  0  0  0 

013  023  033  0  0  0 

^  0  0  0  044  0  0 

0  0  0  0  055  0 

_  0  0  0  0  o  q66 

Since  the  strains  are  defined  in  the  global  axes  (x,  y,  z),  see  Eq.  (6),  it  is  necessary  to 
rewrite  Eq.  (8)  in  that  frame  of  reference.  Specifically,  it  is  found  that 


s  =  R(h)2'  =  R{h)QH-hk-R{h)QQ:'T  0°) 

where  R(‘h  )  denotes  the  6x6  matrix  describing  the  rotation  from  (x',y’,z)  to  (x,  y,  z)  by 
the  ply  angle  <j>*  .  That  is, 


R(h)= 


m 

n2 

0 

mn 

0 

0 


n 


m 


0 
0 

0  1 

-  mn  0 

0  0 

0  0 


-2m  n 
2  mn 
0 

2  2 
m  -n 

0 

0 


0 

0 

0 

0 

m 

n 


0 

0 

0 

0 

-  n 
m 


(11) 


where  m  =  cos§k  and  n  =  sin<j>  k .  For  simplicity  of  notation,  Eq.  (10)  will  be  rewritten  as 

g_-  Q\s- a T\  (12) 


where 

Q  =  R{</>k)QR{~fa )  md  aT  =  [R(-<f>k)]~X  c£T=  R(<t>k)tiT  (13) 
Equations  (l)-(3),  (6),  and  (8)-(13)  provide  a  complete  description  of  the 
displacement,  strain,  and  stress  fields  in  terms  of  the  functions  w0(x,  y,  t),  v0(x,  y,  t), 
w0(x,  y,t),  (p u(x,y,t),  and  cpv(x,  y,().  Then,  to  obtain  the  set  of  governing  equations 
for  these  unknowns,  Hamilton’s  principle  can  be  used,  i.e. 


h 

8  J  (r  -  F  +  Wext)  dt  =  0  (14) 

h 

where  T,  V,  and  Wext  denote  respectively  the  kinetic  and  potential  energies  and  the  work 
done  by  the  external  (acoustic)  loading.  These  quantities  are  readily  obtained  from  the 
displacement,  strain,  and  stress  fields  as 

-  a  b  h/2  ,  \ 

r  =  -JJ  J  p  (m2  +  v2  +  w2]dz  dy  dx  (15) 

2  0  0  -h/2 


(16) 
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where  np  represents  the  number  of  plies  and  zk ,  zk+l  are  the  transverse  coordinates  of 

the  bottom  and  top  of  the  kth  ply.  Finally,  the  external  work  is  defined  by  its  virtual 
counterpart 

a  b 

SJVgxt  =  J  J  p(x,  y,  t )  5 w  dy  dx  (17) 

0  0 

where  p(x,  y,  t )  denotes  the  acoustic  pressure  exerted  on  the  surface  of  the  panel. 

The  application  of  Hamilton’s  principle,  Eq.  (14),  to  the  above  quantities  yields  a 
set  of  five  nonlinear  coupled  partial  differential  equations  for  the  quantities  u0(x,y,t), 
v0(x,y,t),  w0(x,y,t),  <p u(x,y,t),  and  cpv(x,y,/)  the  solution  of  which  represents  a 
serious  challenge  even  in  the  case  of  a  deterministic  loading  p(x,  y,  t).  In  the  presence  of 
random  acoustic  pressure  fluctuations,  the  determination  of  the  statistical  description,  e.g. 
probability  density  functions,  of  the  five  variables  is  beyond  current  capabilities  and  a 
simplification  of  the  model  must  be  achieved. 

2.2  Derivation  of  Simplified  Models 

Several  methods  are  available  for  the  derivation  of  simplified  models  for 
nonlinear  structural  dynamic  problems.  One  standard  approach  is  to  proceed  as  for  linear 
systems  and  express  the  unknown  displacements  fields  in  a  limited  modal-type 
expansion.  In  the  present  context,  this  strategy  would  lead  to  the  approximation 

m 

y)  (18) 

/=1 

where 

7J  =  [«0  (x,  y,  t),  v0  (x,  y,  t),  w0  (x,  y,  t),  <pK  (x,  y,  t),  q>v(x,  y,t)]  (19) 

and  vj/(x,  y)  denotes  a  five-component  vector  of  specified  functions  of  the  spatial 

coordinates  x  and  y.  Finally,  the  time-dependence  of  the  displacement  fields  is  captured 
by  the  unknown  variables  #,(/).  When  the  governing  equations  for  Z  are  linear  or 
weakly  nonlinear,  the  number  of  modes  m  can  be  selected  to  be  substantially  smaller  than 
the  dimension  of  the  vector  Z  and  a  sometime  dramatic  simplification  of  the  problem  is 
accomplished.  A  single-mode  (nr=\)  approximation  is  especially  attractive  as  it  reduces 
the  problem  to  a  single,  nonlinear,  ordinary  differential  equation.  Unfortunately,  in  this 
low  order  approximation  the  components  of  the  vector  Z  are  linearly  dependent  on  each 
other  which  is  not  acceptable  in  the  present  context  since  the  in-plane  displacements 
uq(x,  y,  t )  and  v0(x,  y,  t)  involve  the  transverse  deflections  w0'(x,  y,  t)  only  through  the 

mid-plane  stretching,  a  purely  nonlinear  effect. 

The  derivation  of  a  single-mode  approximation  of  the  present  problem  must  then 
be  accomplished  differently.  Specifically,  it  was  already  argued  in  connection  with  the 
selection  of  the  von  Karman  strains  that  the  in-plane  displacements  are  expected  to  be 
substantially  smaller  than  their  transverse  counterparts.  Thus,  the  contribution  to  the 

panel  kinetic  energy  of  the  terms  w2  and  v2  should  be  small  and  could,  in  first 
approximation,  be  neglected.  Accordingly,  the  displacement  fields  uq(x,  y,t),  v0(x,  y,  t ), 
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(pM  (x,  y,  t ) ,  and  (pv  (x,  y,  t)  appear  explicitly  only  in  the  potential  energy  and  must  then 
be  selected  so  that  8V  =  0.  Note  that  the  potential  energy  does  not  involve  any  time 
derivatives  so  that  its  minimization  with  respect  to  u0(x,y,t),  v0(x,  y,t),  (p u(x,y,t), 
and  cpv(x,  y,  t )  yields  a  set  of  four  nonlinear,  coupled  partial  differential  equations  with 

respect  to  space  for  these  four  displacement  fields.  These  equations  can  be  recognized  as 
integral  versions  of  the  in-plane  equilibrium  condition  for  the  plate  and  can  be  used  to 
express  (at  least  approximately)  uQ(x,y,t),  vQ(x,y,t),  9 u(x,y,t),  and  (pv(x,  y,  t)  in 
terms  of  the  transverse  deflection  w0(x,y,f).  Practically  speaking,  this  procedure 
represents  a  quasi-static  condensation  of  the  plate  equations  and  leads  to  the  single 
unknown  field  w0(x,y,t)  which,  following  previous  arguments  (see  Eq.  (16)-(17»,  can 

be  sought  in  the  form 

wo  (*>  0  =  tf(0  w0  (x,  y)  (20) 

where  w0(x,  y )  is  a  specified  function.  For  maximum  accuracy,  vv0(x,  y)  should  closely 

resemble  the  spatial  distribution  of  the  plate  transverse  deflections  and  thus  can  be 
selected  as  the  panel  buckling  mode  shape.  Mathematically,  this  function  can  be 
expressed  as 

_  /  \  \-i  .  ( WITZ  x\  .  ( Y17T  y\  /Ol'l 

Wo  {x,y)=2^1jamnsm  -  sin  — 7—  (21) 

m  n  v  /  \  / 

where  a  and  b  are  the  dimensions  of  the  plate  in  the  x  andy  directions,  respectively. 

At  this  point  of  the  investigation,  it  is  assumed  that  the  panel  is  simply  supported 
and  thus  the  function  w0  (x,  y)  should  satisfy  zero  deflection  boundary  conditions  on  its 
four  sides,  i.e.  on  x  =  0,  x  =  a,  y  =  0,  and  y  =  b.  These  conditions  are  automatically 
fulfilled  by  the  choice  of  the  sine  functions  independently  of  the  number  of  terms  in  the 
summation.  Thus,  for  simplicity,  it  will  be  assumed  that 


The  above  discussion  demonstrates  that  the  determination  of  a  single-mode 
approximation  of  the  panel  dynamic  behavior  can  be  separated  in  the  following  two  steps: 

(1)  Determination  of  the  functional  dependence  of  u0(x,y,t),  v0(x,  y,t),  cp u(x,y,t), 

and  (pv(x,  y,  t )  on  w0  (x,  y,  t ) 

(2)  Determination  of  the  governing  equation  for  the  transverse  displacement,  i.e.  q(t) . 
Before  addressing  each  of  these  two  aspects  of  the  problem,  the  following  notations  are 
first  introduced: 

Forcesj 

h/2  h/2  hi  2 

NY  =  Joy  dz  My  =  jz  oy  dz  Py  =  J  H(z)  oy  dz  y  =  x,y  (23) 

-h/2  -h/2  -h/2 

h/2  h/2  h/2 

N xy  =  \x xy  dz  M xy  -  Jz  z xy  dz  Pxy  =  J  H{z)  Txy  dz  (24) 

-h/2  -h/2  -h/2 
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N  = 


Nx  Ny  Nxy\ 

hi 2  jjj,  v 

j rf"(2) 


h!  2 


-h/2 


dz 


X z 


dz  Ryz  =  J 


-h/2 


dAyzdz  R 

dz  yz 


Rxz  Ryz 


(25) 

(26) 


Thermal  Loadings 

ht 2~  h/2  ~  h/2 

N‘  =  J  QaT  dz  Ml=  ^zQaT  dz  Pl  =  J  //(z)  Q  aT  dz  (27) 
-h/2  -h/2  -h/2 

where  <2  is  the  3x3  matrix  extracted  from  Q  by  removing  the  3rd,  4th,  and  5th 
rows  and  columns  and  a  is  the  3-component  vector  corresponding  similarly  to  a . 


Strains/Curvatures: 
0 


du0 

1 

+  — 

f 

dv0 

1 

+  - 

(  Bwq  ^ 

dx 

2 

l  dx  J 

dy 

2 

{  dy  ) 

K  - 


a«o_  +  5vo_  + 
dy  dx 

T 


'c/wq  ^ 

dx  J 

dy  J 

n T 


9  9  9 

d  w0  d  wq  2  5  w0 


dxi 


dy1 


dxdy 


$  ~  \fPu  <Pv] 


<!>  = 


d(Pu  dcpv  d(pu  +d(pv 


dx 


dy  dy  dx 


(28) 

(29) 

(30) 


1,  z,z2. 


(dHjzy 

dz 


,zH(z),H2(z),H(z) 


Qdz  .  (31) 


Constitutive  Matrices : 

h/2  f 

(a,  b,  c,  r,  x,w,z)=  I 

-h/2 

Further,  A ,  B ,  C ,  X ,  W ,  and  Z  are  the  3x3  matrices  extracted  from  A,  B,C,  X,  W,  and 
Z,  respectively,  by  removing  the  3rd,  4th,  and  5th  rows  and  columns  and 

i?44  7?45 

_^45  R55 . 

Relations  between  force  and  strain  vectors  can  be  derived  by  multiplying  the 

dH(z ) 


R  = 


(32) 


constitutive  equation,  Eq.  (12),  by  1,  z,  H(z),  and 
panel  thickness.  Specifically,  it  is  found  that 


dz 


and  integrating  through  the 


'N~ 

~A  B  Z  O' 

V" 

n‘ 

M 

B  C  X  0 

K 

m‘ 

P 

Z  X  W  0 

p‘ 

R_ 

0  0  0  R_ 

_  0  _ 

(33) 
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2.2.1  Determination  of  the  functional  dependence  of  Uo(x,y,t),  vq(x,  y, t), 

<p«(*>  y>  and  <Pv(*>  y>  0 on  wo(*>  y>  0 

Proceeding  under  the  assumption  that  the  kinetic  energy  associated  with  the 
displacements/shears  uq (x,  y,  t),  vq (x,  y,  t),  (pu(x,  y,t),  and  <pv(x,  y,t)  can  be  neglected 

in  comparison  with  its  out-of-plane  deflection  component,  the  governing  equations  for 
these  functions  derived  from  Hamilton’s  principle  correspond  in  fact  to  the  minimization 
of  the  potential  energy  alone.  Then,  the  condition  8V  =  0  yields,  after  some  integrations 
by  parts,  the  4  partial  differential  equations 


<5w0  : 

tdN*y_0 
dx  dy 

(34) 

5v0  : 

O 

II 

(35) 

8(pu  : 

o 

11 

1 

g* 

(36) 

3(pv  ■ 

^  ^  Ryz  -  0. 

dx  dy  yz 

(37) 

In  order  to  obtain  closed  form  expressions  for  the  displacements/shears  uq(x,  y,  t), 
v0(x,  y,  t),  q> u{x,y,t),  and  cpv(x,  y,t),  it  is  necessary  to  specify  the  temperature 
distribution  in  the  plate.  Following  Lee  (1993),  it  will  be  assumed  here  that 

TJX  .  7IV  „  m  .  9  7DC  .  y  Tty  1  Z 

T  =  Tq+SvTq  sin—  sin-£  +  Sg  T0  sm2  —  sin2-f--  -  (38) 

u  a  b  *  a  b  4 jh 

where  Tq  denotes  the  average  plate  temperature  and  Sv  and  Sg  represent  measures  of 

the  temperature  gradient  along  and  across  the  plate,  respectively.  Then,  introduce  the 
vectors 
h/2~ 

v  =  j  Q  3  dz 


h/7 

if/  =  jz  Q  adz 

-hi  2  -hi  2 

hi  2 

X=  \H(z)Q5dz 
-hi  2 


1  hn 

6  =  -  \z2  Q  3  dz 

h  1 


-hi  2 


ju  =  -  \zH(z)Q3dz. 
—  h  J 


hi  2 


-hi  2 


(39) 


Under  these  assumptions,  it  is  found  that 


.  7DC  .  7TV 

Tq  +  5V  Tq  sin  sin— — 

a  b 


^  +  \5ST 0 


•  2 m  •  2  1 

sin  — sm - 

a  b  4 


¥ 


(40) 


9 


and 


= 


P*  = 


„  _  _  .  71X  .  7ty 

7q  +  5V  Tq  sin  sin-— 
a  b 


\y  +  SgT0 


.  2  m  ■  iny  1 

sin  — sin - 

a  b  4 
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„  m  .  7zx  . 

T0  +  Tq  sin —  sin— 
a  b 


.  2  ft*  •  2  fty 
sin  — sin  — 

a  b 


(41) 


(42) 


X  +  $g  Tq 

Next,  it  will  be  assumed,  as  is  often  the  case,  that  the  panel  is  composed  of 
symmetrically  placed  layers.  Then,  it  is  readily  seen  that  the  matrices  B,  Z,  B ,  Z ,  and 
vectors  y  and  %  are  identically  zero  and  the  in-plane  displacements  u0(x,  y,  t), 

v0(x,  y,t)  and  out-of-plane  shears  cpM  (x,  y,  t) ,  and  cpv(x,  y,  t)  are  uncoupled  of  each 
other,  see  Eq.  (33)-(37).  The  determination  of  these  two  sets  of  functions  can  thus  be 
accomplished  separately  as  follows. 


(a)  Determination  of  the  in-plane  displacements  uq(x,  y,  t),  vq(x,  y,  t) 

The  determination  of  the  in-plane  displacements  u0(x,y,t),  v0(x,y,t)  is 
accomplished  from  Eq.  (34)-(35)  and  the  reduced  constitutive  relation  Eq.  (33),  i.e. 


N  =  Ae0-Nt,  (43) 

in  either  of  two  ways.  A  first  approach  is  to  introduce  this  expression  for  the  force  vector 
N  in  Eq.  (34)  and  (35)  to  obtain  two  coupled  partial  differential  equations  for  uQ(x,  y,  t), 

vq(x,  y,  t )  (no  time  derivatives).  Another  strategy  is  to  introduce  a  scalar  stress  function 
y)(x,  y,  t )  such  that  (see  Lee  et  al.,  1998  for  details) 


Nx  = 


d  y/ 


Ny  = 


d  y 


o2  a 

N 

*  dxdy 


(44) 


dyA  '  dx z 

so  that  the  equilibrium  equations  (34)  and  (35)  are  automatically  satisfied.  Then,  the 
function  y/(x,  y,  t )  is  selected  to  satisfy  the  compatibility  condition 


d2e°  d2£2 


dy' 


dx‘ 


d2£°3 

dxdy 


f  ~7  ^2 
d  w 


dxdy 


^V,2..A 

dx2 


dy 

dyJ 


(45) 


where  e® ,  e®,  and  £3  are  the  components  of  the  vector  £°,  see  Eq.  (28),  which  are 
related  to  if(x,  y,  t)  through  Eq.  (43)  and  (44).  This  process  results  in  the  4th  order  partial 
differential  equation  for  y/(x,  y,  t)  presented  by  Lee  et  al.  (1998). 

Irrespectively  of  the  approach  selected,  it  is  necessary  to  first  establish  the 
boundary  conditions  to  be  satisfied.  In  view  of  the  simply  supported  nature  of  the  plate,  it 
could  be  expected  that  the  in-plane  motions  of  the  mid-plane  be  restricted  at  the  support, 
i.e. 

u  (0,  y,t)  =  u  (a,  y,t)  =  u  (x,  0,  /)  =  u  (x,  b,  t)  =  v  (0,  y,  t)  =  v  (a,  y,t)  =  v  (x,  0,  t)  =  v  (x,  b,t)  =  0 

(46) 

Of  these  8  boundary  conditions,  it  is  particularly  important  that 

u  (0,  y,  t)  =  u  (a,  y,  /)=  v  (x,  0,  t)  =  v  (x,  b,  t)  =  0  (47) 
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be  satisfied  as  they  correspond  to  the  constraint  on  the  expansion  of  the  plate  and  thus 
play  a  fundamental  role  in  the  buckling  of  the  plate.  The  remaining  four  geometric 
boundary  conditions 

w(x,  0,  t)  =  u(x,b,  /)  =  v(0,y,  t)=v{a,y,t)  =  0  (48) 

represent  a  lack  of  sliding  along  the  supports.  If  necessary  for  simplicity,  they  could  be 
replaced  by  their  corresponding  natural  boundary  conditions 

Txy  {x,  0,  t )  =  TXy  (x,  b,  t )  =  TXy  (0,  y,  t)  =  TXy  {a,  y,t)  =  0  (49) 

on  the  midplane,  z-  0. 

The  determination  of  closed  form  solutions  of  nonhomogenous  partial  differential 
equations,  such  as  those  for  w0(x,  y,  t )  and  v0(x,  y,  t)  or  ft(x,  y,  t),  is  conditional  on  the 
availability  of  both  particular  and  homogenous  solutions.  Considering  first  the  former,  it 
is  seen  that  the  nonhomogenous  character  of  these  equations  originates  with  both  the 
transverse  displacement  w0(x,  y,  t )  and  the  in-plane  temperature  variation  both  of  which 
involve  trigonometric  functions  of  x/a  and  y/b.  On  this  basis,  a  particular  solution  can  be 
sought  as  a  limited  Fourier  series.  On  the  other  hand,  the  determination  of  a  homogenous 
solution  which,  after  superposition  with  its  particular  counterpart,  yields  the  correct 
boundary  conditions  is  in  general  a  daunting  task  unless  the  method  of  separation  of 
variables  can  be  used.  In  this  light,  it  should  be  noted  that 

(i)  the  Fourier  series  particular  solution  for  u0(x,  y,  t)  and  v0(x,  y,  t )  does  not  satisfy 

any  of  the  boundary  conditions  given  by  Eq.  (46)  and/or  (49)  when  and  A2$ 


do  not  vanish. 

(ii)  the  partial  differential  equations  to  be  solved  do  not  admit  a  separation  of 
variables  homogenous  solution  when  A and  A2e  do  not  vanish. 

It  is  then  concluded  that  a  closed  form  solution  for  the  in-plane  displacements  uq(x,  y,  i) 
and  v0(x,  y,t)  will  be  quite  difficult  to  obtain  when  the  two  coefficients  A]6  and  A26 
are  non-zero.  For  this  reason,  the  present  analysis  will  be  limited  to  either  symmetric 
angle-ply  laminates  or  cross-ply  laminates  in  both  of  which  the  constants  A ^  and  A2 g 
always  vanish.  Under  this  restriction,  a  bonafide  solution,  satisfying  Eq.  (47)  and  (49), 
can  be  found  as 


.  2nx 

2  Ttx 

2ny 

- 

(50) 

W01 

sin - h  uq2 

sin 

cos 

a 

a 

b 

.  2rcy 

2kx 

2ny 

(51) 

V01 

sin  +  v02 

b 

cos 

a 

sin 

b 

2  _ 

T  ^  1 

^T0  ,2  t 

*02C 

i  q2 

+  ^•021 

7o)cos 

2  ny 

h 

b 

u 

(52) 

4  tP" 

To 

2  nx 

2ny 

' 

~  2  ^221  5v 

b 2 

cos  • 

a 

cos 

b 
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where 


Ny(x,  y,  /)  =  C50  q2  -  v2  T0  — y-(^200  tf2  +^201  Sv  ?o)cos 

a 


2ftx 


a 


A  ni 


2  7tx  2ny 


^221  Av  7o  cos  -  cos 


9  -’LL  1  -  v  -  u  —  L 

a2  a  b 

/  \  4^2  „  rr,  ‘  2/TX  .  2#^ 

AM  (*>  *  t)  = - —  h2 1  7o  sin -  sin  — 


ab 


a 


(53) 


(54) 


C50  =  — [^22  «2  +  ^12  A2  Qo  =  — y-y  [^12  +  ^1 1  A2 

Sab2  °” u 


Sa2  b2 


hoo  = 


a2  -^li  ^22  ”  A 


12 


32A' 


4l 


^201  =  “ 


a2  (v2  Au  -v,  An) 


\6n‘ 


Au 


^020  - 


b 2  A\  \  A22~  A 


12 


32a‘ 


^22 


^201  = 


b 2  (v'l  ^22  ~v2  An) 


16  ft' 


Vl 


^221  =■ 


16;r 


2  2^r 

W01  =~<1  — 
a 


a22  A\2  ")  Mil  A12 

ZTrvi i 

2  ^11  ,  (A -2  ^12)  ^22 


A  =  ■ 


^22 


^11  ^22  An 
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a2  b2 


+ 


*66 
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\ 


(55) 

(56) 

(57) 

(58) 


^200  ^12 

,2 
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16  A\\  A22  -  A 


1  a2 A22  V\  ~<j2A\2  v2  ~16;t2  A201  A\2 


8  ft  a 


A\  M22  “  A 


Sv  T0 


12 


ft  2  a 
u02  =  77  <?  + 


16#  2  ^4 j  j  ^22  —  ^4 


2 

12  L 


-^(^22  V1  ~A\2  vl)  ~^22\ 
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Mo 
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1  ^020  ^12 


16  aua22-a22 


1  7?2/4u  v2  -  b2A\2  V]  -16;r2  Aq21  A 


8  ftb 


A\\  A22  -An 


1 
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v02  =  777  ^  +  7 - - 
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(b)  Determination  of  the  out-of-plane  shears  qu(x,  y,  t)  and  <pv  (x,  y,t) 

A  set  of  two  partial  differential  equations  for  the  out-of-plane  shears  (pM(x,  y,  t) 
and  cpv  (x,  y,  t)  can  be  derived  from  the  simplified  constitutive  equation,  Eq.  (33),  i.e. 


(63) 


and  Eq.  (36)  and  (37).  An  analysis  of  their  nonhomogenous  terms  motivates  the  search 

for  a  particular  solution  of  the  form 

,  s  r,  .  nx  ny  nx  .  ny 

(pu  {*’  0  =  Use  sm  —  cos  —  +  Ucs  cos  —  sin  — 

0  Cl  0 


'M 

"c 

X 

o' 

K 

Mf 

P 

= 

X 

w 

0 

* 

- 

p‘ 

_R_ 

0 

0 

R 

0 

a 


(  \  . .  >  tv x  ny  T7  nx  .  ny 

(pv  (x,  y,  t)  =  Vsc  sm  —  cos  —  +  Vcs  cos  —  sin  — 


(64) 

(65) 


a  b  a 

t 

Accordingly,  the  vectors  P ,  <p_  ,  and  R  can  be  expressed  as 

nx  .  ny  _  nx 


„  /  >.  _  .  nx  .  ny  n  nx  ny 

P(x,y,t)  =  P  „  sm  —  sin  —  +  PCC  cos — cos  — 
~v  '  a  b  a  b 

>  ,  s'  .  nx  .  ny  '  nx  ny 
<t>  \x,y,t)  =  <j>  ss  sin  —  sin  —  +  <p  Cc  cos — c08” 
“  —  b  ~  a  b 


a 

nx 


a 

nx  .  ny 


„  /  x  „  .  nx  ny  n 

R  (x,y,t)  =  R,r  sin  —  cos—  +  R  cos  —  sm 
_v  )  -sc  b  -cs  b 


(66) 

(67) 

(68) 


Introducing  the  partitioned  vector  Rp  = 
are  equivalent  to 


RT  RT 
—  CS  ~sc 


A SS  E-SS  +Acc  Ecc  Ep  Q 


where 


n  /  a 

0 

0  ‘ 

0 

0 

-n  lb 

0 

0 

n  /  a 

and 

0 

-nib 

0 

Ass  ~ 

0 

0 

nib 

Acc  ~ 

-nl  a 

0 

0 

0 

nib 

0  _ 

0 

0 

-nl  a 

,  it  can  be  seen  that  Eq.  (36)  and  (37) 

(69) 


(70) 


Further,  the  constitutive  relation,  Eq.  (63),  is  equivalent  to  the  equations 

Ess  =  ^  *  *ss  +  ^  £«  ~  Sg  T0  £ 

Ecc  =  Q  X  —CC  +  W  0CC 

where 


(71) 

(72) 


£ss 


2  2 
a2  b 2 


and 


—  SS 


0  0 


2  nl 

"TEb 


t 

Moreover,  from  the  definition  of  (f>  ,  Eq.  (30),  it  is  found  that 


(73) 
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<p_  SS  -  ~  ^SS  —  p  t_  cc  ~  ^ cc  —  p '  0^) 

Then,  combining  Eq.  (69),  (71),  (72),  and  (74)  yields  the  system  of  algebraic  equations 
W  ATss  +  Acc  W  ATCC  +  R2  J Rp  =  q  [ass  X  Ass  +  ^cc  X  £CC]-Sg  T0  Ass  //  (75) 

with 


(76) 


Once  the  coefficients  Usc,  Ucs,  Vsc,  and  Vcs,  stacked  in  the  vector  Rp  have  been 
determined  from  Eq.  (75),  the  particular  solution  for  the  out-of-plane  shears  cpM(x,  y,  t), 
and  (p v(x,y,t)  is  available  from  Eq.  (64)  and  (65).  Further,  the  corresponding  moments 
can  be  evaluated  from  Eq.  (63)  as 

M(x,y,t)  =  Mss  sin  —  sin  —  +  Mcc  cos— -cos—  (77) 

where 

M-ss  =  Q  ^—ss  +  X  ^  ss  ~  3 g  Tq  0  =  q  Ckss  —  X  Ass  Rp  —  Sg  Tq  0  (78) 

M_cc  =  ^  C“icc  +X$_cc=q  Ckcc  -  X  ATcc  R  p  .  (79) 

The  above  derivations  focused  on  the  particular  solution  for  the  out-of-plane 
shears  (p u(x,y,t),  and  q> v(x,y,t).  To  complete  the  solution,  it  remains  to  assess  the 

boundary  conditions  and  to  determine  a  corresponding  homogenous  solution  for  these 
functions.  Considering  first  the  boundary  conditions,  it  is  readily  seen  that  they  do  not 
come  from  the  in-plane  considerations  since  the  contributions  of  cp^(x,  y,  t)  and 
(pv(x,  y,  t )  to  the  displacements  and  strains  vanish  on  the  mid-plane  (z  =  0)  where  the  in¬ 
plane  constraints  are  imposed.  Rather,  the  necessary  boundary  conditions  must  come 
from: 

(i)  the  variational  formulation,  i.e.  5V  =  0 

(ii)  transverse  motion  requirements. 

In  addition  to  Eq.  (36)  and  (37),  the  minimization  of  the  potential  energy  yielded  two 
boundary  terms  which  should  also  vanish.  This  situation  occurs  when: 

either  (pu  is  fixed  or  Px  =  0  and  either  (pv  is  fixed  or  Pxy  =  0  on  x  =  0  and  a  (80) 

and 

either  (pu  is  fixed  or  =  0  and  either  <pv  is  fixed  or  Py  =  0  on  y  =  0  and  b  .  (81) 

Further,  the  simply  supported  nature  of  the  plate  implies  that  the  moments  at  the  supports 
should  vanish,  that  is 

Mx  =0  on  x  =  0  and  a  (82) 


and 

My  =0  on  y  =  0  and  b  .  (83) 

It  is  clearly  seen  from  Eq.  (64)-(66)  and  (77)  that  none  of  the  boundary  conditions  (80)- 
(83)  is  pointwise  satisfied  by  the  particular  solution  derived  above  and  thus  it  should  be 
necessary  to  determine  an  appropriate  homogenous  solution.  The  discussion  conducted  in 
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connection  with  the  homogenous  part  of  the  in-plane  problem  can  be  repeated  here  to 
demonstrate  that  such  a  solution  is  very  unlikely  to  admit  a  closed  form  expression.  A 
first  strategy  to  resolve  this  issue  is  to  restrict  the  current  analysis,  as  in  the  in-plane 
problem,  to  plate  configurations  in  which  boundary  conditions  are  satisfied  exactly  by  the 
particular  solution  presented  above.  This  situation  occurs  if  the  terms  X]6  and  X2(, 

vanish,  as  is  the  case  in  cross-ply  laminates.  Indeed,  it  is  then  found  that  Usc  =  Vcs  =  0  so 
that 

Px  =  0  (pv  =  0  Mx  =0  on  x  =  0  and  a  (84) 

Py=  0  ^„=0  My=  0  ony  =  0  and  b.  (85) 

In  considering  the  boundary  conditions,  it  should  be  noted  that  the  imposition  of  <pv  =  0 
and/or  cpu  =  0  at  the  plate  supports  does  not  imply  that  the  in-plane  displacements  will 

vanish  there  since  the  terms  z  and  z  are  in  general  non-zero.  Thus,  it  is 

dx  dy 

unclear  that  constraining  the  values  of  (pu  and  (pv  at  the  edges  of  the  plate  would  be 
representative  of  a  typical  simple  support  configuration.  Rather,  it  would  appears  that  the 
most  physical  boundary  conditions  for  the  present  problem  would  be 

Px  =  0  Pxy =  0  Mx  =0  on  x  =  0  and  a  (86) 

Py=  0  Pxy  =0  My=  0  ony  =  0  and  b.  (87) 

Further,  note  that  the  above  constraints  are  all  natural  boundary  conditions,  i.e.  they  are 
not  related  to  the  geometry  of  the  problem  but  rather  arise  through  the  variational 
formulation.  Thus,  the  satisfaction  of  these  conditions  is  desired  to  obtain  the  best 
approximation  to  the  exact  solution  but  it  is  not  absolutely  required  as  were  the  in-plane 
constraints,  Eq.  (47).  Finally,  note  that  if  Eq.  (86)  and  (87)  are  not  satisfied  pointwise, 
they  are  matched  in  average  along  the  sides  x  =  0,  a  and  y  -  0,  b.  On  the  basis  of  these 
comments,  it  is  suggested,  as  an  alternative  approach  to  requiring  X jg  and  X2e  to 

vanish,  to  require  that  the  natural  boundary  conditions  be  satisfied  only  on  average  and  to 
let  the  variational  formulation  find  the  best  approximation  possible  to  the  transverse 
response  problem. 

2.2.2  Determination  of  the  governing  equation  for  the  transverse  displacement,  q(t ) 

The  mathematical  developments  presented  in  the  previous  section  have 


demonstrated  that  the  search  for  a  transverse  response  of  the  plate  as 

wo  (*»  0  =  <7 (0  w0  (x,  y)  (88) 

must  be  accompanied  by  corresponding  approximations  of  the  functions  u0(x,y,t ), 
v0(x,  y,  t),  cpw(x,  y,  t),  and  <pv(x,  y,  t )  of  the  form 

wo (x,y,t)  =  q2  «o(x,^)+5v  T0  u0(x,y)+T0  uq (x, y)  (89) 

v0 (x,  y,t)=q2  v0 (x,  y)  +  5 v  T0  v0 (x,  y)  +  T0  v0 (• x ,  y)  (90) 

(*>  y»t)-q  <pM  (x>  y)+5gT< o  9«  (*>  y)  (91a) 

(?v{x,y,t)  =  qyv(x,y)+  Sg  7q  cpv(x,  y)  (91b) 
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where  the  functions  u0(x,y),  ...,  <pv(x,.y)  are  known  from  the  solution  of  the  in-plane 
problem.  Introducing  the  above  expressions  in  Hamilton’s  principle,  Eq.  (14),  results, 
after  some  algebraic  manipulations,  in  the  differential  equation 
f ab 

III 


M  q  +  q  < 


00 


N. 


'  dwr 


\2  f  rsrr.  \2 


dx 


+  N, 
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JJ 

00 


My 


^  32wq  ^ 
dx2 


+  M, 


av2  , 


+  2M 


xy 


^ 82wq  ^ 
dxdy 


dy  dx 


(92) 


-p(t)= o 


where  M  and  p(t)  are  the  modal  mass  and  force,  respectively,  defined  as 
ab  hi 2 

MI  J  p 
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-h/2 
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a  A 


p(0  =  {  J  p(x,  0  ^0  (*>  y)dydx. 


(93) 


(94) 


00 


Note  that  the  expression  for  the  mass  given  in  Eq.  (93)  includes  the  effect  of  the  rotary 

dw0  8w0 

inertia  associated  with  the  pure  bending  rotations  — —  and  — — . 

dx  dy 

Introducing  the  expressions  for  the  forces  N  and  moments  M  given  by  Eq.  (52)- 
(54)  and  (77)  in  Eq.  (92)  yields  the  governing  equation  for  the  transverse  displacements 
in  the  form 

where 
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(96) 
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In  the  above  expressions  (Eq.  (96)  and  (100)),  the  symbols  Mxssq  and  M  yss§  are  the 
first  two  components  of  the  vector  Mss0  which  corresponds  to  the  q  related  terms  of 
Mss  ,see  Eq.  (78).  The  moments  Mxss]  and  MyssX  are  similarly  defined  from  MssX 
which  is  associated  with  the  terms  in  Sg  T0  of  Mss.  Finally,  the  quantities  Mxycc0  and 
Mxycc\  are  obtained  in  a  similar  fashion  from  Mcc ,  Eq.  (79). 


2.2.3  Validation  and  Some  Numerical  Results 

The  mathematical  developments  presented  in  the  two  previous  sections  were 
validated  by  comparison  with  previously  published  results.  In  this  respect,  it  should  first 
be  noted  that  the  present  higher  order  shear  deformation  plate  theory  formulation 
(HSDPT)  naturally  reduces  to  both  the  first  order  shear  deformation  plate  theory 
(FSDPT)  and  the  classical  laminate  plate  theory  (CPT)  under  the  following  assumptions: 

HSDPT  -»  FSDPT:  H(z)  =z  and  R=A  X=W=C  f±  =  6  (101) 

HSDPT  ->  CPT:  H(z)  =  0  and  R=X=W  =  0  //  =  0  .  (102) 


The  present  formulation  was  first  validated  by  comparing  the  expressions  for  the 
parameters  k0,  kx,  k2,  y,  and  p0  obtained  under  the  assumptions  of  Eq.  (102)  for  a 

single  layer  isotropic  plate  with  published  CPT  results  (see  Lee,  1993).  After  some 
algebraic  manipulations,  the  exact  values  were  recovered  as 


kn  = 


n 4 


(a2+62  f 


E 


48 


3  ,3 
a  b 


(103) 


7t4  h 

64 


(104) 

(105) 

(106) 

(107) 


where  a,  E  and  vare  the  coefficient  of  thermal  expansion,  Young’s  modulus  and 
Poisson’s  ratio,  respectively. 

Next,  the  composite  plate  investigated  by  Kavallieratos  (1992)  was  considered 
and  the  present  CPT  results  were  found  to  yield  the  values  of  the  natural  frequency 

(yficolM  Hn=  64.2695  Hz)  and  buckling  temperature  {-k2 lk§  =  9.352°F)  stated  by  the 

author. 

The  validation  of  the  FSDPT  and  HSDPT  formulations  was  accomplished  by 
comparing  the  results  obtained  by  the  present  approach  with  those  published  by  Reddy 
and  Phan  (1985)  for  three  of  their  cases.  Case  I  and  II  relate  to  single  layers  of  a  material 
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that  is  isotropic  (for  case  I)  and  orthotropic  (for  case  II)  while  in  case  III  the  plate  is  a  4 
layer  [079079070°]  composite.  In  all  three  cases,  a  =  b  =  10/j.  Shown  in  Table  1  is  a 
comparison  of  the  natural  frequencies  obtained  by  Reddy  and  Phan  (1985)  and  by  the 
present  formulation.  Note,  in  the  work  of  Reddy  and  Phan,  that  the  inertia  associated  with 
the  in-plane  displacements  u0(x,y,t),  v0(x,y,t)  and  shears  (pM(x,  y,t),  (pv(x,y,t)  is 
included.  Thus,  a  perfect  equality  of  results  should  not  be  expected  but  the  close  matching 
shown  in  Table  2.1  indicates  that  these  inertia  terms  affect  only  slightly  the  transverse 
deflections  as  assumed  here. 

Table  2.1.  Comparison  of  dimensionless  frequencies  obtained  by  Reddy  and  Phan  (1985) 
and  by  the  present  formulation.  Square  plate  with  alh  =  10. 


CPT 

FSDPT 

HSDPT 

Reddy  &  Phan 
(Case  I,  Table  1) 
Present 

0.0955 

0.0955 

0.0930 

0.0934 

0.0931 

0.0930 

Reddy  &  Phan 
(Case  II,  Table  2) 
Present 

0.0493 

0.0492 

0.0474 

0.0476 

0.0474 

0.0473 

Reddy  &  Phan 
(Case  III,  Table  3) 
Present 

18.652 

18.738 

15.083 

15.535 

15.270 

15.054 

Having  validated  the  accuracy  of  the  present  formulation,  it  was  desired  to 
investigate  the  effects  of  plate  thickness  and  shear  deformations  on  the  transverse 
deflections.  It  was  first  observed  that  only  the  coefficients  kQ  and  p$  are  affected  by  the 
choice  of  plate  theory  (CPT,  FSDPT,  or  HSDPT)  as  the  remaining  ones  depend  solely  on 
the  in-plane  stress  field.  Table  2.2  and  2.3  present  the  values  of  k0  and  p0  for  different 
variations  of  the  plate  considered  by  Kavallieratos  (1992)  composed  of  eight  layers 
oriented  at  [07457-45790°]s .  The  geometric  and  material  properties  were  selected  as  a  = 

20  in.,  b  =  8.2  in.,  p  =  1.302 10-4  lb-sec2/in4,  Eu  =18.6  106  psi,  E2 2  =2.0  106  psi„ 
G\2  =  0.8  106  psi,  and  v12  =0.31.  The  values  of  the  shear  moduli  G^  and  G22,  which 
are  required  by  both  the  FSDPT  and  HSDPT,  were  not  specified  by  Kavallieratos  and 
thus  they  were  selected  equal  to  G12.  The  plate  thickness  which  was  selected  by  this 
author  as  h  =  0.0416  in.  (or  blh  =  197.12)  was  varied  here  to  yield  different  ratios  blh, 
specifically  blh  =  197.12,  50,  20,  and  10.  The  different  values  of  k0  and  p0  are 
summarized  in  Table  2.2  and  2.3,  respectively. 
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Table  2.2.  Values  of  the  stiffness  k0  according  to  different  theories 
and  for  different  thicknesses 


b/h 

CPT 

FSDPT 

HSDPT 

197.12 

36.212 

36.206 

36.204 

50 

2218.76 

2212.77 

2210.77 

20 

34668.10 

34091.88 

33903.71 

10 

277344.77 

259844.26 

254498.47 

Table  2.3.  Values  of  the  force  p0  according  to  different  theories 
and  for  different  thicknesses 


b/h 

CPT 

FSDPT 

HSDPT 

197.12 

0.009336 

0.009334 

0.009334 

50 

0.1451 

0.1447 

0.1446 

20 

0.9069 

0.8914 

0.8869 

10 

3.6275 

3.3921 

3.3287 

It  is  seen  from  the  above  tables  that  the  three  theories  yield  similar  values  of  k0  and  pQ 

as  long  as  b/h  remains  larger  than  50  or  so.  As  this  ratio  decreases,  the  difference  between 
the  classical  plate  theory  and  the  shear  deformation  formulations  steadily  increases.  At 
b/h  =  10,  the  classical  plate  theory  is  approximately  10%  away  from  HSDPT  while  the 
two  shear  deformation  formulations  differ  from  each  other  by  only  2%. 

2.3  Panel  Displacement-Stress  Relations 

The  structural  dynamic  modeling  of  composite  panels  accomplished  above  has 
focused  on  the  determination  of  the  differential  equation  for  the  dynamic  response,  i.e. 
the  variable  q(t).  In  addition  to  this  equation,  the  estimation  of  the  fatigue  life  of  these 
panels  requires  also  the  stress-displacement  relations  which  are  obtained  as  follows.  First, 
the  stresses  at  a  certain  depth  z,  a(z),  are  related  to  the  corresponding  strains  s(z) 
through  Eq.  (8).  Next,  the  strains  e(z)  must  be  expressed  in  terms  of  the  single  modal 
displacement  q(t).  This  second  step  is  accomplished  by  relying  on  the  strains,  curvatures 

and  shear  related  terms  s° ,  k,  <p,  and  ^  defined  in  Eq.  (28)-(30).  Specifically,  it  can  be 
shown  that 
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s(z)  =  G{z) 


(108) 


where  the  6x1 1  matrix  G(z )  is 


1 

0 

0 

z 

0 

0 

0 

0 

H(z) 

0 

0 

0 

1 

0 

0 

z 

0 

0 

0 

0 

H(z) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

H'{z) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

H'(z) 

0 

0 

0 

0 

0 

1 

0 

0 

z 

0 

0 

0 

0 

H(z)_ 

(109) 

It  remains  then  to  express  the  four  vectors  e  ,  k  ,  </>,  and  ^  in  terms  of  the 
modal  displacement  q(t).  To  this  end,  note  from  previous  developments,  i.e.  Eq.  (33),  that 

e°  =A~l  N  +  N*  (110) 

in  the  case  of  a  symmetric  composite.  This  relation  provides  a  direct  connection  between 
s°  and  q(t)  since  the  vector  N‘  and  the  matrix  A  are  independent  of  this  quantity,  see 
Eq.  (27)  and  (31),  and  the  components  Nx,  Ny,  and  of  N  are  specified  by  Eq. 

(52)-(54).  Accordingly,  the  vector  £°  can  be  expressed  as  the  sum  of  a  constant  term  and 
a  quadratic  component  in  q,  i.e. 

=£o+§.02<12  ■  OH) 

Turning  next  to  the  curvatures,  k  ,  the  assumed  transverse  displacement 

WQ{x,y)  =  qsm  —  sin—  (112) 

a  b 


(112) 


implies  that 


n  .  nx  .  ny 

——sin —  sin  — 

«2  *  b 

2 

n  .  nx  .  ny 

K-q  —  sin  —  sin  — 

b2  a  b 

2  n  nx  ny 

- cos —  cos - 

a  b  a  b 


(113) 


Finally,  the  shear  related  terms  $  and  <p_  are  readily  evaluated  from  Eq.  (64),  (65),  and 
(67)  where  the  vector  Rp  is  obtained  by  solving  the  system  of  equations  (75).  These 

t 

relations  demonstrate  that  both  (j)  and  <f>  depend  linearly  on  q,  i.e. 
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and 


#  =  $Q+0lq  and  ^  =^0  q ■  (114),  (115) 

Then,  combining  Eq.  (8),  (108),  (111),  and  (113)-(1 15),  it  is  found  that  the  stress 
vector  at  any  point  can  be  expressed  as 


cr  =  cr0  +q_xq  +  a2 


(116) 


2.4  Non-Dimensionalization  of  the  Structural  Dynamics  Equations 

In  order  to  compare  the  effects  of  the  acoustic  excitation  on  different  composite 
and  isotropic  panels,  it  is  convenient  to  first  rewrite  the  governing  equation  for  the 
transverse  displacement  q  and  the  stress-displacement  relations  in  dimensionless  forms. 
For  the  former,  consider  the  differential  equation,  Eq.  (95),  and  introduce  the 
dimensionless  displacement  q  and  time  r  defined  as 

(117) 


-  4 


and 

t  =  0  r  (118) 

where  6  is  an  appropriate  time  constant.  Introducing  these  new  variables  in  Eq.  (95) 
leads  to 


x  e2  ,  a  ^  ,  4  e2  y  /r2  „3  e2p0 

9  +  -*o(i — *)»+— jj— « 


0^ 


where 


s  =  - 


M 

k\  +  ^2 
*0 


2  Mh  2  Mh 


/>(0x) 


To 


(119) 


(120) 


is  the  ratio  of  the  panel  temperature  to  the  buckling  temperature.  Equation  (1 19)  is  of  the 
same  form  (without  the  damping  term)  as  the  one  considered  by  Lee  (1993)  in  connection 
with  isotropic  panels  provided  that 

l2  ^ 


92  2 

Mk 0 


1  + 


a‘ 


Accordingly,  the  time  parameter  6  will  be  selected  as 


0  = 


2 


1  + 


(121) 


(122) 


'X 

Then,  the  nonlinear  term  is  y  q  with 


r  = 


4  d2yh2 
M 


and  the  constant  effect  of  the  temperature  gradients  is 

-  e2^o 

P°  2  Mh' 

Finally,  the  normalized  acoustic  excitation  p(r)  can  be  written  as 


(123) 


(124) 
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(125) 


Mr)= 


2  Mh 


so  that  its  power  spectral  density  is 


S—  (co)  = 

nn\  * 


o 


'pp 


O) 
~0 . 


e 


ppy  '  4  M2h2 


e 


4  M  h 


2  ,2  SPP 


{ co^ 

u 


(126) 


Note  in  the  isotropic  case  that  the  above  expressions  reduce  to  the  relations 
obtained  by  Lee  (1993)  and  used  in  the  previous  reports.  Accordingly,  it  is  desired  that 
the  normalization  of  the  stress  coefficients  also  yields  their  isotropic  forms  which  were 


obtained  by  factoring  the  common  term  E 


n2  h2 


.  Paralleling  this  effort,  let 


°  =  Eeq 


X2  h 2 


(c0+C,«  +  C292) 


(127) 


where  the  coefficient  vectors  C0 ,  Cj ,  and  C2  are  defined  as 

u2 

—  <?<> 


C0=  , 


7t*  Eeq  h‘ 


£i-  .2 


2  b‘ 


n  Eeqh 


^1 


and 


4  y 


E. 


—2 


(128) 

(129) 

(130) 


eq 


To  complete  this  normalization,  it  remains  to  specify  an  “equivalent”  Young’s  modulus 
Eeq .  Relying  on  the  definition  of  k0  for  isotropic  panels,  it  is  suggested  here  that 


^ eq 


3  l3 
a  b 


**  h3(a2+b2J 


(l  -  ^12  ^2l)^0- 


(131) 


2.5  The  Prototypical  Equation 

The  derivations  of  the  previous  sections  have  demonstrated  that  the  structural 
dynamic  response  of  the  panels,  composite  or  isotropic,  is  governed  by  the  same 
equations,  i.e. 

q  +  2i;q  +  (02  (l-s)q  +  yq3  =p0  +  pit)  (132) 

for  the  displacements  and  Eq.  (127)  for  the  stresses.  Accordingly,  most  of  the  ensuing 
discussion  will  focus  on  a  single  example  of  application  without  lack  of  generality. 
Specifically,  the  clamped  isotropic  titanium  plate  already  investigated  by  Vacaitis 
(1994)  will  be  reconsidered  here.  The  dimensions  of  the  panel  were  selected  as  a  =  8.2 
in,  b  =  20.0  in,  and  h  =  0.06  in.  The  acoustic  excitation  was  assumed  to  be  a  white  noise 
random  process  (uncorrelated  fluctuations  over  time  but  fully  correlated  over  space)  of 
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spectral  density  SQ  =  —  105mi°  where  p0  =  2.9 10~9  psi  is  the  reference  pressure 

2n 

and  SPL  denotes  the  sound  pressure  level  in  dB.  The  panel  was  assumed  to  be  subjected 
to  a  uniform  temperature  increase  of  Tq  so  that  the  constant  term  p0  in  Eq.  (132) 
vanishes.  Then,  in  their  dimensionless  forms  (see  Lee,  1993),  the  parameters  of  Eq.  (132) 
and  ( 1 27)  can  be  expressed  as 

“0  +2(J2/3  +  l)  (133) 

s  =  T0  (134) 

Y  =  EPijp2+|3-2+2v  +  i(l-v2[iI()2+p-2)+4(3  +  p-,)"2+(4p+r1)"2+|l+4P"1)' 

3  i  9^8 

(135) 

(136) 
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cos27cx  cos27c;y  cos27tx  cos  4ny  +  cos4rrx  cos27r>1 


(137) 


(138) 


2^  +  P"1)2  (p+4P"]f  4(4p  +  p-1  f 

where  p  =  bla ,  x  =  xla ,  y  =  ylb ,  z  =  zlh.  Further,  the  damping  ratio  £  and  Poisson’s 
ratio  v  were  set  to  0.01  and  0.34,  respectively.  Finally,  in  the  above  dimensionless  form 
the  power  spectral  density  of  the  white  noise  pressure  term  p(t )  is  expressed  as 

(l-v2f  4  Sq 


S—  =144 
PP 


f  b^ 


9  0 


where 


0  = 


(12  1-v2  Jpb4 


(139) 


(140) 


tT  Ehz 

The  numerical  integration  of  the  equation  of  motion,  Eq.  (132)  proceeded  as 
follows.  First,  an  estimate  of  the  time  scale  in  the  fluctuations  of  the  process  q{t)  was 
obtained  by  considering  the  undamped  linear  natural  frequency  co0  and  a  time  step  of 
0.2/coo  was  set  for  the  numerical  integration  of  Eq.  (132).  This  selection  implied  a 
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maximum  (Nyquist)  frequency  (see  Oppenheim  and  Schafer,  1975,  for  a  discussion)  of 
(o^=57ico0  so  that  both  the  acoustic  excitation  p{t)  and  its  corresponding  panel 
response  q{t )  were  defined  in  the  frequency  interval  coef-co^,©^].  Then,  the  constant 
value  of  the  acoustic  loading  in  the  n,h  time  step,  pn ,  was  generated  according  to  its 

specified  power  spectrum  (see  Mignolet,  1993,  for  some  fast  algorithms  to  accomplish 
this  computation).  If  the  excitation  process  is  white  noise  with  constant  power  spectral 

density  Sq  ,  then  the  samples  pn  form  a  sequence  of  independent  random  numbers  with 


common  mean, 


=  0, 


and  variance  E\ 


=  1071(00^0-  Next,  Eq.  (132)  was 


numerically  integrated  with  the  specified  time  step  by  a  Runge-Kutta-Vemer  of  orders 
five  and  six  (IMSL  routine  DIVPRK).  Finally,  the  discrete  values  of  the  response 
deflection  q(t)  were  used  to  produce  the  time  histories  of  the  stresses  according  to  Eq. 


(127). 

The  results  of  this  numerical  integration  are  shown  in  Fig.  2.2-2.13.  A  few 
preliminary  observations  can  be  drawn  from  these  figures.  First,  note  the  clear 
nonlinearity  of  the  displacement-stress  relation  in  the  neighborhood  of  the  bottom 
buckling  state.  Indeed,  although  the  displacement  time  histories  corresponding  to  sound 
pressure  levels  of  119  and  134  dB  are  fairly  symmetric  with  respect  to  this  level,  the 
stresses  are  not,  exhibiting  a  definite  “bottoming  out”.  That  is,  during  the  excursions  of 
the  panel  around  this  position,  the  stress  process  achieves  a  minimum  value  dictated  by 
the  quadratic  relation  (127).  This  peculiarity  is  not  encountered  around  the  top  buckling 
position  since  Eq.  (127)  remains  monotonic  in  this  region. 

A  different  perspective  on  the  behavior  of  the  displacements  can  be  gathered  from 
the  power  spectral  density  plots  of  Fig.  2.8-2.10.  In  particular,  it  is  observed  that  this 
spectrum  exhibits  only  one  peak  at  a  sound  pressure  level  of  1 04dB,  corresponding  to  the 
fluctuations  around  the  buckling  states.  As  the  SPL  increases,  so  does  the  amplitude  of 
the  motions  and  the  response  exhibits  increasingly  nonlinear  features,  as  the  second 
dominant  frequency  shown  in  Fig.  2.9  for  a  SPL  of  119dB.  Note  that  this  second 
frequency  appears  to  be  a  subharmonic  of  order  1/2  of  the  “fundamental”.  Finally,  at  still 
higher  SPL,  see  Fig.  2.10  for  SPL  =  134dB,  the  displacement  process  has  lost  its 
narrowbandedness  and  exhibits  a  single,  rather  wide  peak. 
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Figure  2.4  Time  history  of  the  displacement,  s=  1.8,  SPL  =  134dB 


Figure  2.5  Time  history  of  the  stress,  s  =  1.8,  SPL  =  104dB 
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Figure  2.7  Time  history  of  the  stress,  5=1.8,  SPL  =  134dB 
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Figure  2.1 1  Power  spectral  density  of  the  stress  process,  s  -  1.8,  SPL  -  104dB 
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Density 


Figure  2.12  Power  spectral  density  of  the  stress  process,  s  =  1.8,  SPL  =  1 19dB 


Figure  2.13  Power  spectral  density  of  the  stress  process,  s  =  1.8,  SPL  =  134dB 
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SECTION  3 

EQUIVALENT  LINEARIZATION  TECHNIQUES 

The  estimation  of  the  fatigue  life  of  the  panels  can  be  viewed  as  a  two-step 
process.  Indeed,  it  is  first  necessary  to  evaluate  the  statistical  characteristics  of  the 
dynamic  response,  i.e.  of  the  displacement  q(t )  and  velocity  q{t)  satisfying  the  nonlinear 
stochastic  differential  equation  (132).  In  the  second  stage,  the  damage  generated  by  the 
stresses,  see  Eq.  (127),  is  then  estimated  from  the  known  moments  of  q(t)  and  q{t). 

The  determination  of  the  exact  values  of  the  moments  of  the  displacement  q(t) 
and  the  velocity  q(t)  satisfying  Eq.  (132)  and  corresponding  to  a  white  noise  excitation 
/?(t)  is  easily  accomplished  by  relying  on  the  joint  probability  density  function  of  these 
two  random  variables,  as  given  by  Lutes  and  Sarkani  (1997).  The  availability  of  this 
distribution  is  however  rather  accidental  as  Eq.  (132)  belongs  to  a  limited  class  of 
stochastic  differential  equations  for  which  an  exact  solution  of  the  Fokker-Planck 
equation  can  be  obtained.  For  example,  such  a  closed  form  solution  is  unavailable  if  the 

acoustic  excitation  p(t)  is  a  colored  noise.  Further,  known  expressions  for  the  joint 
distribution  of  the  solution  of  a  system  of  stochastic  differential  equations,  as  would  be 
obtained  for  example  in  a  multiple  mode  analysis  of  the  panel,  are  extraordinarily  limited. 
These  comments  clearly  indicate  that  the  development  of  a  general  strategy  for  the 
estimation  of  the  fatigue  life  of  the  panels  cannot  rely  on  exact  expressions  for  the  joint 
probability  density  function  of  the  displacements  and  velocities.  Rather,  it  is  necessary  to 
rely  on  an  alternate,  general  purpose  methodology  for  the  estimation  of  the  moments  of 
q(t)  and  q(t).  To  this  end,  the  equivalent  linearization  method  (see  Roberts  and  Spanos, 
1 990)  will  be  used. 

According  to  this  methodology,  the  nonlinear  equation  (132)  is  replaced  by  an 
“equivalent  linear”  one  of  the  form 

q  +  2^(o0q  +  keqq  =peq+p{t)  (141) 

where  the  parameters  keq  and  peq  are  selected  so  that  Eq.  (141)  represents  “at  best”  Eq. 
(132).  Specifically,  these  coefficients  will  be  chosen  so  that  the  modeling  error 

£mod  =  E  fcog  (1  -  s)  q  +  y  q3  -  p0  )-  (keq  q  -  peqf  (142) 

where  2s[  ]  denotes  the  operator  of  mathematical  expectation,  is  minimized. 


3.1  Equivalent  Linearization  Strategy  #1 

Proceeding  with  a  differentiation  of  £mocj  with  respect  to  peq  and  keq  yields, 

respectively, 

O>0  (1-*)P<7+Y%3]  -7?0  ~keq  V-q  +  P eq  =  0  (143) 

©0  (l-j)(oJ  +\i2q)+J  Eq4}  -Poflq- keg  P",  \lq=0  (144) 
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where 

a2q=E 


I JLq  and  a 


2  denote  the  mean  and  variance  of  q(t),  i.e.  | xq=E[q]  and 
.  It  is  clearly  seen  from  Eq.  (143)  and  (144)  that  the  evaluation  of  the 


parameters  keq  and  peq  requires  the  knowledge  of  the  yet  unknown  moments  \xq  and 
a2.  This  indetermination  is  resolved  by  replacing  the  exact  values  of  the  mean  \iq  and 

variance  g2  of  the  displacement  q{t)  satisfying  Eq.  (132)  by  those  associated  with  the 

solution  of  the  linear  differential  equation  (141).  Specifically,  following  standard  random 
vibration  arguments  (see  Lutes  and  Sarkani,  1 997),  it  is  found  that 

Pe 

E<7=- 


eq 
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PP 


veq 


2  C  ®0  keq 


(145),  (146) 


Then,  the  solution  of  the  four  coupled  nonlinear  algebraic  equations  (143)-(146)  yields 
the  values  of  \xq,  a2,  peq,  and  keq.  Of  primary  interest  here  are  the  two  moments  and 

thus,  eliminating  peq,  and  keq  from  the  above  equations  yields 
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Note  further  that  the  evaluation  of  the  mean  and  variance  is  accomplished  on  the  basis  of 
the  linear  stochastic  differential  equation  (141)  or  equivalently  under  the  assumption  that 
q{t)  is  a  Gaussian  random  process.  Then,  proceeding  consistently,  the  third  and  fourth 
order  moments  appearing  in  Eq.  (147)  and  (148)  can  be  expressed  as 
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In  solving  Eq.  (147)-(150)  for  the  required  moments  \iq  and  a2,  two  separate  cases  must 
be  considered,  i.e.  \xq  *  0  and  \iq  =0.  Since  the  latter  condition  is  possible  only  when 
=  0 ,  the  analysis  will  focus  separately  on  the  panels  experiencing  a  temperature 
gradient  through  the  thickness  ( Pq  ^  0)  or  a  lack  thereof  (pq  =  0 ). 


3.1.1  Equivalent  Linearization  Results:  Pq^O 

After  some  algebraic  manipulations  (including  a  multiplication  of  Eq.  (147)  and 
(148)  by  \x2),  it  is  found  that  the  mean  value  ]iq  satisfies  the  following  sixth  order 
algebraic  equation 
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Once  the  value  of  pq  has  been  obtained  from  the  above  equation,  the  variance 
evaluated  from  Eq.  (147)  and  (149)  as 
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3y  Eg 


can  be 


(152) 


3.1.2  Equivalent  Linearization  Results:  Pq=0 

When  p0=0,  Eq.  (147)  admits  a  symmetric  solution  of  the  form  1^=0  but  this 
solution  is  not  necessarily  the  only  one  so  that  the  two  separate  sub-cases  p0  =0;\iq  *0 
and  p0  =  0 ;  \iq  =  0  must  be  investigated. 


1.  l.B.  1  Sub-case  # 1 :  pG  =  0 ;  \iq  *  0 

If  the  solution  \xq  =  0  is  not  desired,  the  algebraic  manipulations  carried  out  in 
connection  with  the  derivation  of  Eq.  (151)  can  be  repeated  to  yield 
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so  that 
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Once  the  value  of  \iq  has  been  obtained  from  the  above  equation,  the  variance  aq  can  be 
evaluated  from  Eq.  (152)  as 
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The  above  solutions,  Eq.  (154)  and  (155),  exhibit  some  very  interesting  features.  In 
particular,  in  the  limit  of  a  zero  acoustic  excitation  S--  ->  0 ,  the  mean  value  converges 

toward  the  buckling  states  Q\  and  Q2.  Further,  if  5  <  1,  i.e.  when  the  thermal  effects  are 

not  sufficient  to  buckle  the  plate,  there  is  no  positive  value  of  \iq  in  Eq.  (154)  so  that  no 

solution  of  this  type  exists.  Finally,  the  existence  of  \xq  requires  the  term  inside  the 

square  root  to  be  positive,  or  equivalently  that  (l  -  s)2  -  6  y  — - >0.  This 
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inequality  implies  that  a  nonzero  mean  is  not  possible  when  the  acoustic  level  exceeds  a 
certain  threshold.  Numerical  results  indicate  that  this  threshold  closely  match  the  sound 
pressure  level  (SPL)  at  which  the  snap-throughs  become  very  frequent. 


1.1.B.2  Sub-case  #2:  p0  =  0 ;  \iq  =  0 

On  the  contrary  of  the  mean  value  given  by  Eq.  (154),  the  solution  \xq  -  0  always 


exists  when  p0  =  0  since  Eq.  (147)  and  (149)  are  identically  satisfied.  The  corresponding 
value  of  the  variance  is  then  given  by  Eq.  (148)  and  (150),  that  is 
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3.2  Equivalent  Linearization  Strategy  #2 

It  is  known  (see  Lutes  and  Sarkani,  1997,  for  example)  that  the  exact  probability 
density  function  of  the  displacement  q{t )  is  bimodal  for  all  values  of  the  sound  pressure 
level,  or  equivalently  for  all  values  of  S~~,  provided  that  s>\.  The  equivalent 

linearization  formulation  developed  above  is  consistent  with  this  property  as  long  as  there 
exist  mean  values  satisfying  Eq.  (151)  or  (154).  For  larger  values  of  the  sound  pressure 
level,  however,  the  above  equivalent  linearization  fails  to  accurately  capture  this  property, 
for  example,  for  /?0  =0  and  S—  large  enough  the  only  acceptable  mean  value  \iq  is 

zero  and  the  bimodal  character  is  lost.  To  recover  this  important  property,  it  is  suggested 
here  to  proceed  with  a  different  equivalent  linearization  formulation  in  which  the  mean  is 
imposed  to  be  different  from  zero  and  the  standard  deviation  is  obtained  to  minimize  the 
modeling  error  of  Eq.  (142).  Physical  arguments  motivate  the  selection  of  the  imposed 
mean  values  to  be  the  buckled  states,  i.e.  \iq  -  Q\  or  \x.q  -  Qj  >  where 


ooq  {l-s)Qi  +  yQj  =0  Qi*  0  /  —  1,2. 


Accordingly,  Eq.  (145)  requires  that 
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so  that  the  modeling  error,  Eq.  (142),  can  be  rewritten  as 
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Differentiating  the  above  expression  with  respect  to  the  remaining  parameter,  keq ,  yields 


the  equation 
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Then,  using  Eq.  (146),  (149),  and  (150),  it  can  be  shown  that  the  above  relation  reduces  to 
the  following  equation  for  the  variance  of  q, 
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the  solution  of  which  is 
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Note  that  this  solution  always  exists  for  5  >  1  as  expected  from  the  bimodal  character  of 
the  exact  probability  density  function. 


3.3  Numerical  Results 

The  prototypical  problem  analyzed  in  the  previous  section  was  reconsidered  here 
to  assess  the  validity  and  accuracy  of  the  equivalent  linearization  strategies.  For 
simplicity,  it  was  again  assumed  that  the  effects  of  the  gradient  through  the  thickness  are 

negligible  so  that  pG  =  0 . 

Then,  shown  in  Fig.  3.1  are  the  standard  deviations  of  the  displacement  q  obtained 
by  Monte  Carlo  simulation  of  the  fully  nonlinear  equation  Eq.  (132)  and  by  the 
equivalent  linearization  strategies  #1  and  #2  for  5  =  1.8  as  a  function  of  the  sound 
pressure  level  ( SPL )  in  dB.  The  behavior  of  the  Monte  Carlo  results  is  particularly 
informative,  at  very  low  SPL  the  panel  vibrates  with  a  very  low  amplitude  around  one  of 
the  buckled  states.  Since  these  two  positions  are  equally  probable  and  equally  distant 
from  the  undeformed  position  q  =  0,  the  overall  mean  is  zero  and  the  corresponding 
standard  deviation  should  be  very  close  to  Q\ ,  as  confirmed  by  Fig.  3. 1  (for  5  =  1.8,  Q\  = 
0.448).  As  the  sound  pressure  level  increases  however,  two  different  phenomena  occur. 
First,  the  level  of  vibration  around  the  buckled  states  increases  slightly  and 
unsymmetrically,  i.e.  the  vibrations  are  larger  toward  the  undeformed  position  because  of 
the  decreased  local  stiffness  exhibited  by  the  restoring  forces.  This  lack  of  symmetry  of 
the  response  induces  a  decrease  of  the  mean  displacement  of  the  motions  around  each  of 
the  buckled  states,  i.e.  pj  <Q\ .  The  second  effect  is  the  appearance  of  snap-throughs, 
although  very  infrequent  at  first,  that  populate  the  region  in  between  the  two  buckled 
states.  Both  of  these  factors  imply  the  decrease  of  the  standard  deviation  seen  in  Fig.  3.1. 
This  process  continues  until  the  response  becomes  dominated  by  the  snap-throughs  in 
which  case  an  increase  in  the  excitation  level  mainly  induces  an  increase  in  the  level  of 
response,  i.e.  the  maximum  displacement  away  from  the  undeformed  position. 
Correspondingly,  the  standard  deviation  then  starts  increasing  as  a  function  of  the  sound 
pressure  level  as  shown  in  Fig.  3.1.  These  observations  would  tend  to  confirm  the  earlier 
conjecture  that  the  motion  at  low  SPL  is  governed  by  the  fluctuation  processes  only  while 
at  high  excitation  levels  it  is  the  snap-throughs  that  dictate  the  panel  response.  More 
importantly,  it  would  appear  that  the  region  in  which  both  of  these  processes  are 
important  is  quite  narrow  and  is  located  near  the  minimum  of  the  standard  deviation. 
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Table  3.1.  Maximum  value  of  the  sound  pressure  level  for  which  nonzero  mean 
equivalent  linearization  solutions  exist  as  a  function  of  temperature  s. 


5 

1.05 

1.8 

3 

5 

SPL  (dB) 

97  1 

120 

129" 

135 

Considering  next  the  equivalent  linearization  results,  note  first  that  when  nonzero 
mean  solutions  exists,  the  corresponding  overall  standard  deviation  can  be  computed  as 

a<7  +  +^2  +  a*).  (163) 

Next,  on  the  basis  of  physical  arguments,  the  equivalent  linearization  results 
corresponding  to  the  first  approach  were  obtained  for  “low”  SPL,  i.e.  for  sound  pressure 
levels  lower  than  the  threshold  values  given  in  Table  3.1,  by  the  combination  of  the 
models  associated  with  nonzero  mean  values  \xq  closest  to  Q\  and  Q2.  However,  when 

the  excitation  was  large  enough  that  no  real  solution  of  Eq.  (154)  could  be  obtained,  the 
zero  mean  model  characterized  by  the  variance  of  Eq.  (156)  was  used.  The  results  of  Fig. 
3.1  clearly  demonstrate  that  the  combination  of  the  two  nonzero  mean  solutions  yields  an 
extremely  accurate  estimate  of  the  overall  standard  deviation  but  that  its  zero  mean 
counterpart  (for  SPL>  120dB)  can  severely  underestimate  this  moment  although  the 
accuracy  appears  to  be  improving  as  the  sound  pressure  level  increases  past  the 
fluctuation  to  snap-through  transition  region.  Figure  3.1  also  demonstrates  that  the  second 
equivalent  linearization  strategy  is  not  as  reliable  as  the  first  one  except  for 
SPL  >  120  dB.  In  fact,  it  exhibits  the  wrong  trend  in  the  low  sound  pressure  level  regime, 
the  estimated  standard  deviation  increases  as  a  function  of  the  SPL.  This  undesirable 
property  is  associated  with  the  fixed  value  of  the  means  p,-  =Qt  which  does  not  allow 
the  good  modeling  of  the  unsymmetry  of  the  local  stiffness  discussed  above.  On  the  basis 
of  these  results,  it  is  suggested  to  use  the  equivalent  linearization  strategy  #1  when 
nonzero  means  exist  and  to  rely  the  second  strategy  when  such  solutions  are  not  possible. 
The  above  trends  occur  consistently  through  the  investigated  range  of  temperatures 
se[l,  5]. 

It  should  further  be  noted  from  Fig.  3.1  that  the  average  of  the  standard  deviation 
estimates  corresponding  to  the  equivalent  linearization  strategies  1  and  2  approximates  in 
fact  remarkably  well  its  exact  counterpart  for  SPL>  120.  A  similar  accuracy  has  also 
been  observed  at  other  values  of  the  temperature  s  indicating  that  this  averaging  could 
consistently  be  used  to  refine  the  equivalent  linearization  estimates. 

A  different  perspective  on  the  reliability  of  the  equivalent  linearization  techniques 
can  be  obtained  by  comparing  the  produced  probability  density  functions  with  their  exact 
counterparts.  This  comparison  is  presented  in  Fig.  3.2-3.7.  Note  first  that  the  improved 
accuracy  of  the  standard  deviation  estimate  obtained  from  the  equivalent  linearization  #1 
at  “low”  SPL  is  obtained  by  allowing  a  mean  shift  of  the  corresponding  distribution.  This 
effect  is  small  but  apparent,  see  Fig.  3.2  for  SPL  =  114dB  and  s  —  1.8.  This  shift 
increases  with  the  sound  pressure  level  and  provokes  an  undesirable  error  in  the 
prediction  of  the  location  of  the  peaks  of  the  probability  density  function,  see  Fig.  3.3  for 
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SPL  =  119dB.  The  equivalent  linearization  strategy  #2,  on  the  contrary,  maintains  the 
location  of  the  peaks  at  the  buckling  states  as  shown  in  Fig.  3.4  but  overestimates  the 
height  of  the  peak.  Note  also  that  both  methods  fail  to  capture  the  probability  of  crossing 
the  origin,  i.e.  of  the  occurrence  of  snap-throughs. 

At  higher  SPL  still,  the  equivalent  linearization  strategy  #1  fails  to  produce  a  non¬ 
zero  mean  solution  resulting  in  a  substantial  overshoot  of  the  probability  of  crossing  the 
origin,  see  Fig.  3.5.  Although  the  second  equivalent  linearization  approach  appears  to 
capture  quite  well  the  exact  probability  density  function  for  this  condition,  as  shown  in 
Fig.  3.6,  it  also  will  produce  single  peak  distribution  as  the  sound  pressure  level  continues 
to  increase  or  equivalently  as  the  temperature  decreases,  as  can  be  noted  from  Fig.  3.7. 
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Figure  3.1 .  Standard  deviations  of  the  response  as  functions  of  the  sound  pressure  level 
obtained  by  simulation  and  by  the  equivalent  linearization  strategies  #1  and  #2  for  5  =  1 .8 . 


Figure  3.2  Probability  density  functions  of  the  displacement,  exact  (— )  and  estimated 
according  to  the  equivalent  linearization  #1  (+++),  s  =  1.8,  SPL  =  1 14dB. 
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Figure  3.4  Probability  density  functions  of  the  displacement,  exact  (— )  and  estimated 
according  to  the  equivalent  linearization  #2  (□□□),  s  =  1.8,  SPL  =  1 19dB. 
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Figure  3.5  Probability  density  functions  of  the  displacement,  exact  (--- )  and  estimated 
according  to  the  equivalent  linearization  #1  ( —  ),  s  =  1 .8,  SPL  =  134dB. 


Figure  3.6  Probability  density  functions  of  the  displacement,  exact  (--- )  and  estimated 
according  to  the  equivalent  linearization  #2:  around  top  buckling  state  (  —  )  and 
combined  (□□□),  s  =  1.8,  SPL  =  134dB. 


40 


SECTION  4 

FATIGUE  DAMAGE  PREDICTION 

The  second  predictive  aspect  of  this  investigation  focused  on  the  estimation  of  the 
accumulated  damage  by  using  the  equivalent  linearization  results.  That  is,  it  was  desired 
to  duplicate  at  best  rainflow  results  (Downing  and  Socie,  1982)  obtained  from  the 
numerically  evaluated  displacement/stresses  time  histories  from  the  mean  and  variances 
given  by  Eq.  (151),  (152),  (154)-(156),  and  (162).  To  this  end,  three  approaches  of 
increasing  complexity  were  investigated,  all  of  which  rely  in  some  fashion  on  the 
Rayleigh  approximation  (see  Lutes  and  Sarkani,  1997)  and  on  an  expected 
narrowbandedness  of  the  response  processes  (displacement,  velocity,  stresses,  etc.).  In 
assessing  the  properties  of  each  of  these  approximations,  it  should  first  be  noted  that  the 
displacement  q  is  not  a  Gaussian  process  because  of  the  cubic  nonlinearity  of  Eq.  (132) 
and  that  the  displacement-stress  relationship  is  also  nonlinear  as  seen  in  Eq.  (127).  Thus, 
even  if  the  displacements  were  Gaussian,  the  stresses  would  not.  In  this  light,  the  three 
formulations  to  be  presented  in  sections  4.2-4.4  rely  on  the  following  assumptions: 

(1)  both  stresses  and  displacements  are  Gaussian  (“standard”  Rayleigh 
formulation) 

(2)  the  displacements  are  Gaussian  but  the  stresses  are  not,  the  nonlinear 
displacement-stress  relationship  describes  the  distribution  of  stresses 
(nonlinear  displacement-stress  formulation) 

(3)  the  displacements  are  specified  as  the  sum  of  Gaussian-type  processes 
describing  the  specificities  of  the  fluctuations  around  the  buckled  states  and 
the  snap-throughs.  In  this  detailed  model,  the  nonlinear  displacement-stress 
relationship  describes  the  distribution  of  stresses  (phenomenological 
formulation). 

Before  these  three  distinct  approaches  are  described,  however,  an  exact  formula 
for  the  accumulated  damage  will  be  presented  that  was  found  quite  useful  in  interpreting 
the  numerical  results. 


4,1  Damage  Accumulated:  An  Exact  Formula 

In  assessing  the  damage  accumulation  in  the  panel,  it  is  first  assumed  that  the 
material  is  characterized  by  the  S-N  curve 


Nf=KS;m  (164) 

where  N  f  is  the  number  of  cycles  to  failure  when  the  stress  range  is  Sr ,  and  K  and  m 

are  material  constants.  Further,  adopting  a  linear  damage  accumulation  rule,  the  total 
damage  in  the  panel  after  a  time  T can  be  estimated  as 

.  ,  Wc(7>«)  1  ,  nATfm) 
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where  nc[Tfm)  denotes  the  total  number  of  half-cycles  in  the  time  interval  t  e  [0  ,Tfm 


and  Sr  k  represents  the  stress  range  in  the  klh  half-cycle.  Note  that  this  quantity  can  also 
be  represented  as 

Sr,t  =  \s(lk ) - d  =  1 S(t)a| -  (‘Js(r)|A  (166) 

where  tk_\  and  tk  denote  the  beginning  and  ending  times  of  the  stress  range  Srk  . 
Further,  the  last  equality  in  the  above  relation  holds  because  of  the  constant  sign  of  the 
stress  velocity  S(t)  in  the  interval  te[tk_ \,tk].  Combining  Eq.  (165)  and  (166)  in  the 


special  case  m  =  1  yields 
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Moreover,  relying  on  the  stationarity  of  the  stress  velocity  process,  it  is  found  that  the 


expected  damage  accumulated  in  the  time  interval  t  e  [0,  T j-m 


can  be  expressed  as 


(is*) 

A  final  simplification  of  the  above  relation  can  be  obtained  by  noting  from  Eq.  (127)  that 

S  =  q  (Cj  +  2  C2  q)  (169) 

where  the  velocity  q  is  known  to  be  Gaussian  with  mean  zero  and  standard  deviation 


a, 


and  to  be  independent  of  the  displacement  q  (see  Lutes  and  Sarkani, 


1997).  Then, 


e[\s\}=  E[\q\]E[\Cl+2C2  q\ ]=  E[\C\  +2C2  4 


(170) 


and  the  expected  damage  is  given  by 

£[b(7>„)]=-= L-  of  E[\CI+2C2  ?|]  Tfm.  (171) 

A  generalization  of  the  above  relation  for  m  *  1  is  unfortunately  unavailable,  even 
for  integer  values  of  m  .  Indeed,  in  these  cases,  an  expression  for  the  expected  damage  in 
terms  of  the  m  point  correlation  of  the  stress  velocity  process  can  be  derived  by 
proceeding  as  in  Eq.  (166)-(168).  However,  the  lack  of  existence  of  a  closed  form 
solution  for  the  transition  probability  density  function  of  the  displacement  q  and  velocity 
q  prevents  the  derivation  of  a  manageable  expression  for  the  expected  damage.  One  must 
then  resort  to  either  rainflow  simulations  or  to  the  use  of  approximate  relations  as  derived 


in  the  next  sections. 
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4.2  Damage  Accumulated;  “Standard”  Rayleigh  Formulation 

The  simplest  approximate  expression  for  the  accumulated  damage  used  in  the 
present  investigation  is  traditionally  referred  to  as  the  Rayleigh  approximation  (see  Lutes 
and  Sarkani,  1997).  It  is  derived  by  assuming  that 

(i)  the  stress  process  S(t)  is  Gaussian 

(ii)  the  stress  process  S(t)  is  narrowband  (required  unless  m=  1). 

Under  these  two  conditions,  it  can  be  shown  (see  Lutes  and  Sarkani,  1997)  that 
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(172) 


where  <js  and  cr^  are  the  standard  deviations  of  the  stress  and  stress  velocity  processes. 
Further,  E[AD]  denotes  the  expected  damage  accumulated  over  one  cycle  of  the  stress 
process  which  can  be  estimated  as 
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It  is  known  that  the  Rayleigh  approximation  is  exact  when  m  =  1  under  weaker 
conditions  that  the  two  stated  above,  (i)  and  (ii).  In  fact,  when  m  =  1,  it  is  only  necessary 
that  the  stress  velocity  process  be  Gaussian.  This  property  can  be  confirmed  by  noting 
that  Eq.  (174)  reduces  to  Eq.  (171)  when  m  =  1  and  C2  =  0 . 

The  determination  of  the  standard  deviations  ct^  and  can  be  accomplished 
from  Eq.  (127)  and  (169).  Specifically,  relying  on  the  Gaussian  character  of  q  and  q ,  it  is' 
found  that 

=C\  °2q+C2  [4V2q  +2°2q)Gl  +4C1  C2  Vq  °q  (175) 


and 

CT^  =G4  1C12  +4C2  ^4  +a^)+4Cl  C2  Vq\  •  (176^ 

It  should  be  noted  that  this  formulation  is  simple  but  it  is  also  inconsistent:  it  assumes  that 
both  the  stresses  and  displacements  are  Gaussian  and,  at  the  same  time,  that  they  are 
nonlinearly  related  through  Eq.  (127)  and  (169). 


4.3  Damage  Accumulated:  Nonlinear  Displacement-Stress  Formulation 

The  removal  of  the  inconsistency  of  the  Rayleigh  formulation  described  above 
can  be  accomplished  by  formulating  the  entire  problem  in  terms  of  the  displacement  q 
which  here  will  be  assumed  to  be  a  Gaussian  process  (or  a  combination  thereof).  The 
approach  presented  below  follows  the  non-Gaussian  correction  scheme  introduced  by 
Lutes  and  Sarkani  (1997).  For  clarity  of  the  presentation,  assume  momentarily  that  the 
displacement-stress  relationship  S  =  g{q)  is  monotonic  (this  is  clearly  not  always  the 
case  for  Eq.  (127))  and  that  the  response  processes  are  narrowband.  Then,  it  can  be 
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argued  that  the  positive  ( S+ )  and  negative  ( S~ )  peak  stresses  of  a  given  stress  range 
correspond  to  peak  values  of  the  displacement  q.  Further,  under  the  narrowbandedness 
assumption,  the  peak  deviations  from  the  mean  displacement  follow  a  Rayleigh 
distribution  and  the  positive  and  negative  maximum  deviations  of  a  given  cycle  are 
approximately  equal.  That  is, 

+  and  S~  =  g{pq  -aq  u)  (177) 

where  u  is  a  standard  Rayleigh  random  variable.  Accordingly,  the  expected  damage 
accumulated  over  a  cycle  is 

S+  -5“  +aq  u)-s(vq  -°q  U)T  1 "  exp(-  W2  ll\du 

-I  *oL  J 

(178) 

Equation  (178)  is  not  directly  applicable  here  since  Eq.  (127)  is  not  monotonic. 
Indeed,  there  exists  an  extremum  (minimum  or  maximum)  of  the  stress  for 

r2 

C  L, 

q-a0- - —  and  the  corresponding  stress  is  Se  =Cq  -  — — — .  On  physical  grounds, 

2  C2  4  C2 

it  is  expected  here  that  the  membrane  stresses  yield  a  stretching  effect  so  that  C2  >  0  and 
the  extremum  is  always  a  minimum.  Thus,  for  all  values  of  u>ue=  qe  -  \xq  /  a  g,  the 

stress  process  does  not  undergo  a  single  cycle  of  magnitude  S+  -S~  but  rather  two 

separate  cycles  of  respective  amplitudes  S+  -Se  and  Se-S  =  S  —Se  . 
Accordingly,  the  expected  damage  accumulated  over  the  cycle  of  displacement  is 


E[AD]  =  ^-  §\g(\Lq+Oqu)-g(»q-<*q«)\m  U  CXP(- U2  /  2)  du 

K  oL  J 

+  T  J  \s(vq+c<iu)-se\m  +\g{vq-°qu)-se\m  u  exp(-  u2  12}  du  (179) 
K  J I 

where 

S  =  «fe)=c0+C,?  +  C2?2.  (180) 

Finally,  relying  on  the  expected  narrowband  character  of  the  displacement  process 
and  the  assumed  Gaussian  distribution  of  the  displacement  and  velocity  processes,  the 
period  can  be  estimated  through  the  expected  rate  of  upcrossing  of  the  mean  value  as 

T  = - —  so  that  the  expected  damage  accumulated  in  the  time  interval  /  6  [0,  Tjjn  j  is 
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For  example,  for  m  =  1,  Eq.  (179)  becomes 
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(183a-c) 


(184) 


As  a  check,  consider  the  case  C2  =0  for  which  qe=<&  so  that  Eq.  (182)  reduces  to 


£[A£>]=Cj  fin—  and  4°(  "F= ~pr  Oq  Tfm  as  expected  from  Eq.  (171). 

K  ^2%  A 

The  case  C\  =  0 ,  =  0  is  also  very  important  as  it  corresponds  to  the  most 

strongly  nonlinear  displacement-stress  relationship  and  thus  represents  a  good  test  of  this 
second  damage  accumulation  formulation.  Accordingly,  it  is  found  that  qe  =  0  and 


Se=  0  so  that 


E[AD]  =  -C2o 


and  E[D{Tfl J=- 


2  C, 
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4.4  Damage  Accumulated:  Phenomenological  Formulation 

The  nonlinear  displacement-stress  formulation  presented  in  the  last  section 
accounted  fully  and  exactly  for  the  nonlinear  q  to  S  transformation  so  that  improvements 
over  this  strategy  require  the  consideration  of  non-Gaussian  displacements.  In  keeping 
with  the  stated  goals  of  this  investigation  to  obtain  an  estimation  strategy  of  the  damage 
accumulated  that  can  be  extended  to  multi-modes,  it  is  proposed  here  to  introduce  a 
formulation  based  on  a  “combination”  of  Gaussian  processes  that  is  valid  for  panels 
statically  buckled,  i.e.  for  s>l.  Specifically,  in  accord  with  the  non-zero  mean 
equivalent  linearization  strategies  developed  above,  it  will  be  assumed  that  the  motion 
around  each  buckled  configuration  follows  a  Gaussian  distribution  so  that  the  probability 
density  function  of  q  is  the  weighted  sum  of  two  Gaussian  distributions  of  means  exactly 
or  approximately  equal  to  the  buckled  states.  That  is, 
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where  the  probabilities  (weights)  q\  and  q2  are  selected  so  that 
q\+q2=\  and  ©q  (l-s)E[q]  +  yE 
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i.e.,  so  that  q  has  a  total  probability  of  1  and  that  Eq.  (132)  is  satisfied  at  least  on  average. 
In  particular  when  p0  =  0,  it  is  directly  found  that  q\  -q2  =1/2  as  expected.  Note  in 
Eq.  (186)  that  the  means  pj ,  p2  and  standard  deviations  cj ,  a2  can  be  estimated  from 
the  equivalent  linearization  #1  (provided  solutions  pj ,  p2  to  Eq.  (151)  exist)  or  #2. 

Once  the  probability  density  function  of  q  ,  Eq.  (186),  is  known,  it  is  necessary  to 
evaluate  the  damage  accumulated.  To  this  end,  it  is  suggested  here  to  investigate 
separately  the  fluctuations  around  each  buckled  configurations  and  the  snap-throughs 
from  one  such  position  to  the  other.  Both  fluctuation  and  jump  processes  will  be 
considered  to  be  narrowband  so  that  their  corresponding  damages  can  accurately  be 
estimated  by  Rayleigh’s  formula, 

E{D{Tfm)\  =  nE[kD]  (188) 

where  n  is  the  expected  number  of  cycles  occurring  in  the  time  interval  of  length  Tfin .  To 
exemplify  the  determination  of  the  damage  per  cycle,  E[ED\,  and  the  number  of  cycles  n, 
consider  a  half-cycle  of  motion  with  q(t k-\)^  k)  and  note  that  there  are  three  distinct 

possibilities  for  the  displacement  time  history,  i.e. 

(a)  fluctuation  around  the  top  buckled  state 

(b)  snap-through  from  top  to  bottom  buckled  state 

(c)  fluctuation  around  the  bottom  buckled  state 

as  depicted  in  Fig.  4.1.  To  evaluate  the  likelihood  of  each  of  these  three  trajectories,  it  is 
first  necessary  to  estimate  the  probabilities  p\  and  p2  that  the  peak  displacement 
q[tk_ i )  is  around  the  top  and  bottom  buckled  states,  respectively.  To  this  end,  assume  for 

simplicity  that  ~pQ  =  0  and  introduce  the  probability  of  snap-through  p  so  that  the 

probabilities  of  the  trajectories  (a)-(c)  are 

pa=p\(\-p)  Pb=P\P  ^  Pc =  Pi-  (189a-c) 

Accordingly,  the  probability  that  the  valley  of  the  half-cycle  be  located  near  the  bottom 
buckled  state  is  p ^  +  pc  which  should  also  equal  to  p\  since  the  displacement  process  q 

is  symmetric  when  pQ  =  0 .  It  is  then  required  that 

Pb  +Pc  =P\  P  +  Pl  =Pl-  (19°) 

Since  the  total  probability  for  the  peak  location  must  equal  1,  one  also  has  the  condition 

P\+P2=1  (191) 

and  the  solution  of  Eq.  (190)  and  (191)  yields 

p  ,=  —  and  (192a,b) 

2-p  2,-p 

In  fact,  the  peak  to  valley  transition  process  has  been  represented  as  a  Markov 
chain  with  a  transition  probability  p  which  need  now  be  evaluated.  This  determination 
requires  a  physical  characterization  of  the  occurrence  of  snap-throughs.  To  this  end,  note 
that  snap-throughs  are  associated  with  the  existence  of  a  softening  region,  see  Fig.  4.2,  in 
which  the  local  stiffness  is  negative.  The  domain  of  attraction  around  one  of  the  buckling 
states  could  be  defined  as  the  region  in  which  the  potential  is  less  than  or  equal  to  the 
value  at  q  =  0.  The  values  Qyp  and  Q2T  bounding  this  region  are  then  such  that 
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(193) 
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from  which  it  is  found  that 
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(194) 


Then,  it  is  suggested  that  a  snap-through  occurs,  i.e.  trajectory  (a)  turns  into  (b),  when  the 
extreme  displacement  corresponding  to  the  valley  of  the  half-cycle,  q(tb  ),  falls  outside  of 
the  domain  of  attraction,  i.e. 

r- 

{Q\t  -Mi)2 
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where  the  last  equality  holds  in  view  of  the  expected  Rayleigh  distribution  of  the  peak 
displacement  values  q{tb ). 

The  combination  of  Eq.  (192)  and  (195)  yields  the  probabilities  associated  with 
the  three  possible  trajectories  (a)-(c).  Further,  the  corresponding  damage  contributions 
2?[AZ)fl],  F,[ADbl  and  £[ADC]  can  be  evaluated  as  in  Eq.  (179),  (180),  and  (182)  so  that 

the  total  damage  can  be  evaluated  as 

E[D{Tfm)\=na  E[M)a]+nb  E[ADb]  +  nc  E[ADC}  (196) 

where  na ,  nb,  and  nc  are  the  number  of  cycles  corresponding  to  each  type  of  trajectory 
(a)-(c).  These  values  can  be  expected  to  be  proportional  to  the  probability  of  occurrence 
of  their  corresponding  trajectory,  Eq.  (189),  that  is 

na  =pi  (\-p)N  nb  =  p\  p  N  and  nc  =  p2  N  (197) 

where  N  denotes  the  total  number  of  cycles.  In  turn,  this  quantity  can  be  estimated  from 
the  total  time  T pm  as 
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where  the  upcrossing  rates  va,  vb,  and  vc  are  estimates  of  the  frequencies  of  the 
narrowband  trajectories  (a)-(c).  The  Gaussian  character  of  the  motions  around  the 
buckled  states  suggests  that 
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and 
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are  the  standard  deviations  of  the  velocities  associated  with  the 

fluctuations  around  the  top  and  bottom  buckled  states  and  can  be  determined  from  either 
equivalent  linearization  #1  (provided  solutions  pj,  p2  to  Eq.  (151)  exist)  or  #2.  In  fact, 


where  a^  and 


they  are  both  equal  to  the  exact  expression  a  •  = . 
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The  determination  of  proceeds  similarly  but  in  connection  with  the  overall 
probability  density  function  of  the  displacements,  i.e.  Eq.  (186).  Specifically,  the 
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upcrossing  rate  of  the  undeformed  panel  configuration,  <7  =  0,  is  given  in  terms  of  the 
joint  probability  density  function  of  q  and  q  as  (see  Lutes  and  Sarkani,  1997) 

Vb  =  v0  =  jV  Pqq  (°> V)  dv  ■  ^200) 

0 

Assuming  that  the  displacement  q  and  velocity  q  are  statistically  independent  and  that 
the  latter  random  variable  is  Gaussian,  it  is  found  that 

V4=V»  =p<i{0)'jh  <201) 

where  pq(0)  is  given  by  Eq.  (186)  for  q  =  0 .  Note  the  assumptions  stated  above  are  not 

particularly  restrictive  since  they  are  known  to  hold  for  the  processes  satisfying  either  the 
exact  equation  of  motion,  Eq.  (132),  or  any  of  its  equivalent  linearization  approximations. 

Combining  Eq.  (197)-(199)  and  (201)  yields  estimates  of  the  number  of  cycles 
na,  nb,  and  nc  corresponding  to  each  type  of  trajectory  (a)-(c)  and  completes  the 
damage  accumulation  formulation. 


4.5  Numerical  Results 

To  establish  baseline  values  the  accumulated  damages,  the  numerical  integration 
of  the  equation  of  motion  (132)  was  performed  as  discussed  in  section  2.  From  the  time 
histories  of  the  displacement,  the  corresponding  realizations  of  the  stress(es)  were 
obtained  through  the  quadratic  relation,  Eq.  (127).  Then,  the  rainflow  cycle  counting 
strategy  of  (Downing  and  Socie,  1982)  was  used  to  identify  stress  ranges  and  evaluate  the 
accumulated  damage  over  a  fixed  interval  of  Tfm=  32,000  time  steps.  In  keeping  with  the 


dimensionless  formulation  of  the  equation  of  motion,  the  damage  was  normalized 

\m 


according  to  D  =  K 


KEeq  n2  h2  J 


D .  Finally,  the  values  of  the  damage  presented  in 


these  figures  are  in  fact  averages  over  100  realizations  of  stress  time  histories  to  reduce 
the  variability  of  the  damage  estimates.  These  100  simulations  were  divided  into  two 
groups  of  50  each  with  the  excitation  records  in  the  second  set  equal  in  magnitude  but  of 
opposite  sign  to  their  counterparts  in  the  first  set.  This  procedure  was  required  to  ensure 
that  exactly  half  of  the  simulations  at  low  sound  pressure  levels  would  lead  to  fluctuations 
around  the  upper  buckling  position  with  the  remaining  half  around  the  lower  one. 

Shown  in  Fig.  4.3-4.10,  are  the  normalized  damages  computed  on  the  basis  of  the 
stress  ax  at  the  middle  of  either  the  top  surface  of  the  panel  (Z  =  0.5)  or  its  neutral  plane 
(Z=  0).  Physically,  it  can  be  seen  that  this  stress  component  is  dominated,  when  Z=  0.5, 
by  the  linear  bending  terms,  C\  q ,  which  is  completely  absent  when  Z  =  0.  Accordingly, 


it  can  be  expected  that  this  latter  situation  provides  a  worst  case  scenario  from  the 
standpoint  of  nonlinearity  of  the  displacement-stress  relations. 

Turning  now  to  the  prediction  approaches,  the  equivalent  linearization  strategy  #1 
was  used  to  provide  the  estimate  of  the  standard  deviation  of  the  displacement  process 
required  by  the  “standard”  Rayleigh  and  nonlinear  displacement-stress  formulations.  On 
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the  contrary,  the  phenomenological  approach  relied  on  the  results  of  the  second 
equivalent  linearization  strategy,  consistently  with  the  basis  of  this  numerical  technique. 
An  analysis  of  these  results  demonstrates  that  both  the  Rayleigh  and  nonlinear 
displacement-stress  formulations  yield  reliable  estimates  of  the  accumulated  damage  at 
“low”  sound  pressure  levels,  i.e.  when  snap-throughs  are  rare  events.  These  estimates 
were  obtained  by  averaging  the  damages  corresponding  to  each  of  the  two  nonzero  mean 
models.  In  effect,  this  procedure  fully  accounts  for  the  fluctuation  processes  but  neglects 
the  damage  associated  with  the  infrequent  snap-throughs. 

As  the  SPL  approaches  the  threshold  above  which  the  equivalent  linearization 
strategy  #1  fails  to  yield  a  nonzero  mean,  snap-throughs  are  becoming  more  frequent  and 
both  the  Rayleigh  and  nonlinear  displacement-stress  formulations  underestimate  the 
actual  damage.  In  this  condition,  the  phenomenological  formulation  provides  very 
reliable  approximation  of  the  accumulated  damage  provided  that  the  value  of  the 
probability  density  function  pq{  0),  governing  the  rate  of  snap-throughs,  is  accurate. 

Finally,  at  high  sound  pressure  levels,  i.e.  higher  than  the  above  threshold,  it  was 
found  that  the  Rayleigh  approach  typically  overestimated  the  exact  accumulated  damage 
and  that  the  nonlinear  displacement-stress  formulation  underestimated  this  value. 
However,  the  average  of  these  two  estimates  (shown  in  Fig.  4.3-4.10  as  “Final”)  has 
consistently  been  found  to  be  an  accurate  approximation  of  the  damage. 

In  view  of  the  above  results  and  comments,  it  is  suggested  that  the  nonlinear 
displacement-stress  formulation  be  used  until  the  sound  pressure  level  approaches  the 
threshold  at  which  the  equivalent  linearization  method  #1  stops  yielding  nonzero  mean 
solutions.  Then,  in  a  small  region  below  this  value,  the  phenomenological  formulation 
should  be  used  with  an  accurate  estimate  of  pq( 0).  Finally,  above  the  threshold,  the 

approximate  accumulated  damage  should  be  obtained  as  the  average  of  the  estimates 
provided  by  the  Rayleigh  and  nonlinear  displacement-stress  formulations. 

The  above  discussions  have  focused  on  the  prediction  of  the  accumulated  damage 
associated  with  a  specific  stress  component.  In  general  however,  the  panel  experiences  a 
multiaxial  state  of  stress  and  it  is  necessary  to  define  an  appropriate  equivalent  stress. 
Adopting  a  Tresca-type  failure  criterion,  it  is  suggested  here  to  compute  the  damages 
associated  with  the  stress  components  ox,  o y ,  and  a  x  -  a y .  The  accumulated  damage 

will  then  be  selected  as  the  largest  of  these  three  values.  This  process  is  demonstrated  in 
Fig.  4.13-4.15  in  connection  to  a  composite  panel  similar  to  the  one  introduced  by 
Kavallieratos  (1992)  and  discussed  in  Section  2  but  with  dimensions  a  =  b  =  8  in.  and 
different  lay-ups.  The  damages  shown  in  Fig.  4.13-4.15  correspond  to  the  point  at  the 
middle  of  the  panel  on  its  top  surface  (Z=0.5). 
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Figure  4.1  Probability  density  function  of  the  displacement  showing 
the  three  peak  to  valley  trajectories  (a)-(c). 


SOFTENING 


Figure  4.2  Force  vs.  displacement  curve  showing  the  softening  region. 
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Figure  4.3  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 
for  m  =  1, 5=  1.8,  Z  =  0.5. 
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Figure  4.4  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 

for  m  =  1,  s-  1.8,  Z- 0. 
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Figure  4.5  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 
for  m  =  2,  s=  1.8 ,  Z  =  0.5. 


Sound  Pressure  Level  (SPL,  dB) 


♦ 

Rainflow 

■ 

Rayleigh 

▲ 

N.L. 

-x- 

Final 

Figure  4.6  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 

for  m  =  2,  s-  1 .8,  Z  =  0. 
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Figure  4.7  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 
for  m  =  5,  s=  1.8,  Z=  0.5. 
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Figure  4.8  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 
for  m  =  3,s=  0.5,  Z  — 0.5. 
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Figure  4.9  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 
for  m  =  3,  s=  1 .05,  Z  =  0.5. 
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Figure  4.10  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 
for  m  =  2>,s-  1.8 ,Z=  0.5. 
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4.1 1  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 

for  m  =  3,  s=  3,  Z  =  0.5. 
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Figure  4.12  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  and  by  the  various  approximate  methods 

for  m  =  3,  s=  5,  Z-  0.5. 
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Figure  4.13  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  for  the  lay-up  [90  45  -45  0]s  with  m  =  3,s=  1.8 ,Z=  0.5 
and  for  the  stresses  ax  (Sx),  cy  (Sy),  and  ox  -  oy  (Sx  -  Sy). 
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Figure  4.14  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  for  the  lay-up  [-45  90  45  0]s  with  m-3,s=  1.8 ,Z  -  0.5 
and  for  the  stresses  ox  (Sx),  cr^,  (Sy),  and  ox  -  a y  (Sx  -  Sy). 
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Figure  4.15  Estimates  of  the  nomalized  damage  as  functions  of  the  sound  pressure  level 
obtained  by  rainflow  analysis  for  the  lay-up  [45  -45  90  0]s  with  m  =  3,s=  1.8 ,  Z  =  0.5 
and  for  the  stresses  ax  (Sx),  Gy  (Sy),  and  gx  -oy  (Sx  -  Sy). 
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SECTION  5 
SUMMARY 

The  focus  of  this  investigation  has  been  on  the  prediction  of  the  fatigue 
life/damage  of  composite  panels  subjected  to  an  extreme  environment,  i.e.  to  both  high 
thermal  effects  (temperature  and  temperature  gradients)  and  a  strong  acoustic  transverse 
loading.  To  achieve  this  goal,  it  was  necessary  to  accomplish  the  following  tasks: 

(i)  develop  an  appropriate  structural  dynamic  modeling  of  the  panel 

(ii)  derive  reliable  approximate  expressions  for  the  statistics  of  the  panel 
response  to  the  random  acoustic  loading 

(iii)  formulate  a  prediction  strategy  of  the  accumulated  damage  in  terms  of  the 
obtained  statistics  of  the  response 

All  three  problems  were  successfully  addressed.  First,  a  large  displacement  -  small  strains 
structural  dynamic  model  of  the  composite  panel  was  accomplished  by  relying  on  the  von 
Karman  strain  expressions.  Further,  the  formulation  naturally  accounts  for  uniform 
temperature  effects  as  well  as  the  presence  of  in-plane  and  transverse  temperature 
gradients.  Finally,  a  higher-order  displacement  field  was  adopted  to  accurately  capture 
the  shear  effects.  Consistently  with  the  proof  of  concept  aspect  of  this  Phase  I  effort,  a 
simplified,  one-mode  approximation  of  the  response  was  derived  for  a  simply  supported 
panel.  Non-dimensionalization  of  the  resulting  equation  of  motion  revealed  that  the 
response  of  different  types  of  panels  should  exhibit  similar  characteristics  and  a 
prototypical  panel  was  considered.  The  analysis  of  its  response  revealed  in  particular  that 
the  panel  response  is  fairly  narrowband  at  low  sound  pressure  levels  when  the  panel 
vibrates  around  its  buckled  states.  However,  as  the  excitation  level  is  increased,  the 
nonlinear  effects  increase  as  well  resulting  in  the  appearance  of  a  subharmonic  (of  order 
1/2)  component  of  the  response.  Finally,  at  very  high  sound  pressure  levels,  the  power 
spectral  density  of  the  panel  response  exhibits  a  single  broad  peak.  The  nonlinearity  of 
the  displacement-stress  was  emphasized  and  a  “bottoming  out”  effect  in  the  stresses  was 
observed. 

The  determination  of  reliable  approximations  of  the  statistics  of  the  panel 
response  to  the  random  acoustic  excitation  was  accomplished  by  relying  on  two  separate 
equivalent  linearization  strategies.  The  first  of  these  two  methods  seeks  to  approximate  at 
best  the  panel  response  by  a  Gaussian  process  of  unknown  mean  and  variance.  At  low 
SPL,  nonzero  means  are  indeed  found  that  correspond  to  the  fluctuations  around  the 
buckled  states.  Above  a  critical  SPL,  however,  this  approach  only  yields  a  zero  mean 
approximation  of  the  “continuous”  snap-through  of  the  panel.  The  second  equivalent 
linearization  method  relies  on  a  mean  equal  to  the  buckled  states  and  a  representation  of 
the  panel  response  is  obtained  as  the  sum  of  zero  mean  fluctuations  around  these 
positions.  A  reliable  approximation  of  the  variance  of  the  response  was  obtained,  at  low 
SPL,  by  the  first  of  these  approaches  alone  but,  at  high  SPL,  a  reliable  estimate  of  this 
quantity  is  only  obtained  by  averaging  the  predictions  of  the  two  equivalent  linearization 
strategies.  Even  though  accurate  approximations  of  this  moment  were  obtained,  it  was 
shown  that  the  exact  probability  density  function  of  the  panel  response  can  accurately  be 
matched  by  sums  of  Gaussian  distributions  only  at  low  SPL. 

The  prediction  of  the  accumulated  damage  in  the  panel  was  achieved  by  a 
combination  of  three  methods:  a  “standard”  Rayleigh  approximation,  a  more  detailed 
formulation  taking  into  account  the  nonlinearity  of  the  displacement-stress  relation,  and 
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finally  a  phenomenological  modeling  of  the  fluctuations  around  the  buckled  states  and 
the  snap-through  process.  When  the  snap-throughs  occur  only  exceptionally,  the 
nonlinear  displacement-stress  formulation  yields  an  excellent  approximation  of  the  exact 
damage  as  estimated  by  an  extensive  simulation/rainflow  analysis.  However,  as  the  sound 
pressure  level  is  increased  and  snap-throughs  start  occurring,  this  approach 
underestimates  the  damage  and  the  phenomenological  formulation  ought  to  be  used. 
When  the  SPL  is  high  enough,  i.e.  above  the  threshold  at  which  the  equivalent 
linearization  strategy  #1  fails  to  yield  nonzero  means,  the  average  of  the  damages 
predicted  by  the  Rayleigh  and  the  nonlinear  displacement-stress  formulations  was  shown 
to  be  quite  accurate  over  a  broad  range  of  sound  pressure  levels  and  S-N  curve  exponents. 

The  combination  of  these  three  different  research  efforts  provides  a  solid 
foundation  for  the  prediction  of  the  fatigue  life  of  composite  and  isotropic  panels  of 
various  shapes,  support  conditions  and  in  a  wide  array  of  extreme  environments. 
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